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EXECUTIVE  SUMMARY 


OBJECTIVE 

The  objectives  of  this  study  were  to  develop  a  basic  understanding  of  the  release  rate  and 
subsequent  fate  and  transport  of  munition  constituents  (MCs)  in  water  and  in  sediment.  The  ability  to 
characterize,  assess,  and  predict  potential  MC  source  loading  and  distribution  has  significant 
implications  for  Department  of  Defense  (DoD)  and  Department  of  Navy  (DoN)  range  sustainability 
initiatives.  DoD  will  gain  critical  information  for  making  scientifically  defensible  risk  management 
decisions  about  underwater  ordnance  leave-in-place  (LIP)  mitigation  and  blow-in-place  (BIP)  vs. 
removal  options.  In  addition  to  explosive  blast  (safety)  considerations,  future  regulatory  emphasis 
will  likely  require  an  assessment  of  potential  underwater  ordnance  contamination  and  mitigation 
efforts  that  could  include  water  and  sediment  quality  issues. 

TECHNICAL  APPROACH 

As  a  baseline  scenario,  this  study  directly  addresses  the  amount  of  MC  introduced  into  the 
environment  from  the  case  of  a  single  breached  munition.  The  information  determined  for  MC 
release  from  a  single  breached  shell  was  incorporated  into  the  calibrated  hydrodynamic/transport 
model  to  predict  fate  and  transport  of  the  MC  and/or  degradation  products.  To  address  the  specific 
release  function  and  its  predictive  ability,  we  conducted  analytical,  empirical,  and  numerical 
modeling  studies  of  MC  source  release  from  breached  shells  under  various  hydrodynamic  and  shell 
integrity  conditions.  Release  rates  were  developed  and  validated  for  a  single  breached  shell  in  three 
nominal  release  scenarios:  a  shell  on  the  surface  of  bottom  sediment,  a  shell  from  a  low-order 
detonation  above  the  bottom  in  the  water  column,  and  a  shell  buried  in  sediment.  The  semi- 
analytical  release  rate  function  was  validated  by  results  from  both  the  empirical  study  and  the 
numerical  modeling  study,  using  the  FLUENT  model.  The  calibrated  model,  TRIM2D,  and 
Environmental  Protection  Agency’s  (EPA)  WASP?  model  were  both  used  for  simulating  fate  and 
transport  of  trinitrotoluene  explosive;  2-methyl- 1,3, 5-trinitrobenzene  (TNT)  released  by  design  from 
a  breach  shell  in  San  Diego  Bay,  California,  and  Elizabeth  River,  Virginia. 

RESULTS 

The  release  function  developed  semi-analytically  assumes  that  the  release  through  a  breach  in  a 
munition  casing  can  be  determined  by  the  following  five  key  parameters:  1)  the  growth  of  the  breach 
hole  (as  the  radius  of  the  hole,  (b);  2)  the  radius  of  the  cavity  formed  due  to  loss  of  mass  released 
from  inside  the  shell,  (R );  3)  the  dissolution  rate  (p)  from  solid  to  aqueous  phases  of  the  MC  inside 
the  shell  casing);  4)  the  outside  ambient  current  (U)  to  which  the  casing  hole  is  exposed;  and  5)  the 
hydrodynamic  dispersion  coefficient,  (D).  This  function  dictates  the  release  rate  be  characteristically 
governed  by  two  Reynolds  Numbers:  the  current-based  Reynolds  Number  (Ub/D),  and  the 
dissolution-based  Reynolds  Number  (pR/D).  Release  rate  is  governed  by  current  (U)  and  dissolution 
rate  (p)  for  the  scenarios  (pR  /D  »n  Ub/D)  and  (pR  /D  «  Ub/D),  respectively.  Empirical  studies 
were  conducted  for  release  scenarios  both  in  the  water  column  and  in  sediment.  For  a  current- 
controlled  scenario,  release  rates  are  proportional  to  the  ambient  current  speed,  and  therefore,  release 
rates  in  the  water  column  are  much  greater  than  in  sediment.  However,  total  time  for  a  complete  MC 
release  from  a  test  shell  (L  =  24”,  r  =  8”,  8-kg  MC  weight)  also  depends  on  the  breach  hole  size  (e.g., 
crack  area),  and  the  dissolution  rate,  which  can  result  in  the  total  release  time  varying  from  a  few 
days  to  several  hundred  years. 
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The  fate  and  transport  modeling  study  shows  that  predicted  TNT  concentrations  are  on  the  order  of 
nanograms  per  liter  (ng/L)  in  the  vicinity  of  the  release  locations  for  a  single  test  shell.  This  provides 
the  basis  for  an  order  of  magnitude  estimate  for  worst-case  scenarios.  For  example,  to  generate  TNT 
concentrations  to  the  micrograms  per  liter  (pg/L)  level  in  the  water  column,  it  would  take  one 
thousand  shells  of  the  same  shell  integrity  to  be  co-located  at  a  single  site. 

BENEFITS 

A  basic  understanding  of  processes  and  governing  factors  for  the  release  rate  and  fate/transport  of 
MC  in  marine  environment  was  developed.  Predictive  modeling  capabilities  were  also  developed  for 
these  processes,  which  will  help  DoD  scientifically  address  the  MC  issues  mentioned  above.  With 
the  data  developed  from  this  study,  DoD  will  be  better  equipped  to  make  technically  defensible 
managerial  decisions  for  sites  with  underwater  ordnance. 
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1.  INTRODUCTION 


This  effort  focused  on  providing  the  information  necessary  for  assessing  fate  and  transport 
processes  associated  with  munition  constituents  (MCs)  release  from  a  single  breached  shell  into  a 
seawater  environment.  The  objective  is  to  provide  a  predictive  modeling  capability  of  MC  fate  and 
transport  associated  with  an  unexploded  breached  projectile.  The  results  are  being  promulgated  to 
assist  Department  of  Defense  (DoD)  and  regulatory  agencies  with  a  science-based  predictive  tool  to 
estimate  the  amount  of  MC  released  to  the  environment  from  underwater  munitions  over  time. 

Our  initial  approach  (Shell  Model  Phases  I  &  II)  focused  on  a  non-empirical  analytical  and 
modeling  scenario,  wherein  a  single  breached  round  is  hypothetically  introduced  into  the  marine 
environment  and  subsequent  MC  release,  fate,  and  transport  are  predicted.  To  understand  the  release 
behavior,  we  have  developed  a  set  of  specific  semi-analytical  formulae  that  describe  the  release  rate 
of  munitions  constituents  under  various  conditions  of  shell  integrity  and  hydrodynamic  situations. 
Specifically,  the  MC  release  rate  can  be  explicitly  expressed  as  a  function  of  the  following  five 
parameters:  ambient  current  speed,  hydrodynamic  difflisivity  coefficient,  size  of  the  breach  hole, 
cavity  volume  inside  the  shell,  and  dissolution  rate  of  MC  from  the  solid  to  aqueous  phase  inside  the 
shell.  Predictions  based  on  this  semi-analytical  release  rate  function  have  been  partially  verified  by 
results  from  a  generic  fluid  dynamic  model,  FLUENT.  Generally,  this  release  rate  function  is 
considered  applicable  for  release  of  any  type  of  chemical  from  breached  shells,  with  breach  openings 
ranging  from  tiny  holes  to  fully  exposed  (low-order  detonation),  under  various  aquatic  environments. 
However,  to  fully  verify  the  semi-analytical  solution,  further  modeling  and  empirical  studies  would 
help  increase  confidence  in  the  underlying  assumptions.  Subsequent  to  model  verification,  a  more 
complete  model  can  be  developed  to  include  integration  of  empirical  MC  release  rates,  shell 
corrosion  rates,  and  MC  rates  of  decomposition.  Use  of  these  empirical  process  descriptors  along 
with  analytically/empirically  verified  release  rate  formulae  will  adequately  describe  and  predict 
source  release  and  fate  and  transport  of  MC  from  a  breached  shell  in  water  and  sediment 
environments,  and  can  ultimately  be  used  for  evaluating  risk  to  the  environment  in  which  the 
breached  shell  resides.  Extrapolation  to  release  by  multiple  shells  is  possible  given  additional 
information  relating  to  the  amount  of  MC  associated  with  the  number,  type,  and  integrity  of 
munitions  in  a  prescribed  area. 

As  an  initial  approach  to  coupling  the  release  model  described  above  to  fate  and  transport,  we  have 
evaluated  two  additional  fate  and  transport  models,  TRIM2D  and  WASP7.  Details  of  the  release  rate 
formulae  development,  analysis  and  evaluation,  and  partial  verification  using  FLUENT,  empirical 
process  descriptors,  and  preliminary  fate  and  transport  modeling  results  using  TRIM2D  and  WASP7 
are  described.  The  WASP7  model  simulates  trinitrotoluene  explosive.  The  2-methyl- 1,3,5- 
trinitrobenzene  (TNT)  concentrations  in  both  water  column  and  sediment  with  box-like  resolution, 
whereas  TRIM2D  provides  simulated  TNT  concentrations  in  water  column  at  much  finer  resolutions. 
Both  models  provide  meaningful  complementary  results  that  are  mutually  consistent,  demonstrating 
their  utility  for  developing  an  understanding  of  MC  transport  in  aquatic  environments. 

Two  estuarine  scenarios  are  evaluated:  San  Diego  Bay,  California,  and  the  Elizabeth  River, 
Virginia.  Due  to  unique  characteristics  of  these  two  estuarine  systems,  the  simulated  transport  and 
flushing  patterns  of  TNT  in  these  two  Navy  harbors  exhibit  significance  differences.  Generally, 
flushing  in  San  Diego  Bay  is  driven  by  tides,  whereas  flushing  in  Elizabeth  River  is  driven  by  both 
tidal  forcing  and  riverine  flows.  It  was  observed  in  the  model  that  differences  in  TSS  settling 
velocities  and  in  the  associated  sediment  bed  is  much  more  pronounced  for  Elizabeth  River  than  for 
San  Diego  Bay,  resulting  in  higher  predicted  TNT  concentrations  in  the  Elizabeth  River  sediments. 
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The  Shell  Model  Phase  III  study  focused  on  verification  of  the  release  rate  function  developed  in 
the  initial  effort  associated  with  MCs’  release  from  a  single  breached  shell  into  a  shallow  seawater 
environment.  The  objective  was  to  provide  a  validated  predictive  modeling  capability  of  MC  release 
rate  associated  with  a  single  unexploded  breached  projectile.  Validated  results  are  required  to  assist 
DoD  and  regulatory  agencies  with  a  science-based  predictive  tool  that  is  technically  defensible  when 
using  the  models  to  estimate  the  amount  of  MC  released  to  the  environment  from  underwater 
munitions  over  time. 

Phase  III  addressed  the  results  of  two  primary  tasks,  the  numerical  modeling  study  and  the 
corresponding  empirical  study  and  four  subtasks.  The  general  fluid  dynamic  model,  FLUENT,  was 
used  to  simulate  the  release  rate  of  MC  leached  from  a  single  breached  shell  under  various 
hydrodynamic  (e.g.,  currents)  and  shell  integrity  (e.g.,  sizes  of  breach  hole  and  cavity)  conditions. 
The  simulated  scenarios  were  designed  and  focused  on  the  following  four  sub-tasks;  1)  the  effect  of 
shapes  of  the  opening,  2)  effect  of  orientation  of  the  crack  opening,  3)  internal  circulation,  and  4) 
probabilistic  modeling.  The  first  three  scenarios  were  simulated  for  a  single  shell,  and  the  last 
scenario  addresses  both  a  single  shell  and  a  group  of  shells  associated  with  probabilistic  functions  for 
currents  and  hole  sizes.  Results  from  these  studies  are  consistent  with  the  analytical  solutions 
previously  developed  in  the  Phase  I  Study. 

The  empirical  study  was  conducted  using  the  flume  tank  at  the  University  of  Wisconsin- 
Milwaukee.  A  surrogate  dye  in  solution,  used  to  mimic  MC  solution,  with  mixing  in  the  cavity  of  a 
breached  shell  by  the  ambient  currents,  were  measured  by  both  PIV  (Particle  Imagery  Velocimetry 
and  LIF  (Laser  Induced  Fluorescence)  instruments.  From  this  empirical  study,  it  was  found,  among 
other  things,  that  release  rate  is  governed  by  the  Reynolds  number  (Ub/D),  the  release  velocity  (Uf), 
and  the  shear  velocity  (u*).  Generally,  the  release  velocity  can  be  correlated  with  the  shear  velocity. 
Correlation  of  the  release  velocity  to  the  shell  internal  cavity  size  is  stronger  than  that  to  the  breach 
hole  size,  particularly  for  a  large  cavity.  All  of  the  empirical  results  compare  well  with  the  numerical 
modeling  results,  further  verifying  the  release  rate  function. 

For  the  initial  three  phases  of  study  described  generally  above,  we  evaluated  MC  release  from  a 
breached  shell  sitting  on  the  sediment  bed,  subject  to  hydrodynamic  mixing  and  transport  by  ambient 
current.  This  resulted  in  the  development  and  validation  of  the  following  MC  release  rate  function, 
for  which  the  derivation  will  be  described  later,  incorporating  five  variables,  breached  hole  size, 
internal  cavity  radius,  ambient  current,  mixing  coefficient,  and  dissolution  rate; 

_  ItiDC^ 

2D  D 
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In  Phase  IV,  we  analyzed  the  order  of  magnitude  of  the  pore-water  advective  transport  and  the 
relative  magnitude  between  low  pore-water  advection  and  the  difflisivity  (mixing)  coefficient.  MC 
release  for  breached  shells  buried  in  sediment  is  dictated  by  the  same  release  rate  function,  except 
that  ambient  current  is  replaced  by  pore-water  flow,  which  generally  is  much  slower  than  the  current 
in  the  water  column.  The  relative  magnitude  of  these  two  quantities  determines  how  the  MC  is 
released  from  a  buried  breached  shell  and  how  MC  in  the  plume  is  transported  in  this  situation. 

The  analysis  of  pore  water  and  its  effects  on  MC  release  and  transport  from  breached  shells  buried 
in  sediment  was  validated  by  both  an  empirical  study  and  two  modeling  studies.  The  empirical  study 
was  conducted  to  measure  release  of  dye  (surrogate  for  MC)  from  a  breached  shell  buried  under  the 
sediment  in  the  “sediment  flume  box.”  Pore-water  advection  in  sediment  within  the  sediment  flume 


2 


box  was  generated  by  maintaining  a  constant  hydraulic  head  difference  at  the  two  ends  of  the  box. 
Pore-water  advection  speed  was  measured  and  compared  with  a  pore-water  advection  model  based 
on  Darcy’s  Law.  The  rate  of  dye  release  from  the  breached  shell  surrogate  buried  in  sediment  was 
recorded  and  then  compared  to  the  theoretical  solution,  which  follows  a  first  order  decay  rate.  The 
pore-water  model  was  further  used  to  characterize  the  relative  length  scales  of  horizontal  and  vertical 
domains,  within  which  pressure-induced  pore-water  advection  occurs. 

As  demonstrated  previously,  the  relative  magnitude  of  pore-water  advection  versus  mixing 
coefficient  (or  diffusion  coefficient  in  sediment)  can  significantly  influence  MC  release  and  transport 
behaviors  in  sediment.  A  model  has  been  constructed  to  examine  the  effects  of  these  two  parameters 
on  MC  plume  transport  in  sediment. 
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2.  BACKGROUND 


The  ability  to  characterize,  assess,  and  predict  potential  munition  constituent  (MC)  source  loading 
and  distribution  has  significant  implications  for  Department  of  Defense  (DoD)  and  Department  of 
Navy  (DoN)  range  sustainability  initiatives.  DoD  would  gain  critical  information  for  making 
scientifically  defensible  risk  management  decisions  about  underwater  ordnance  leave-in-place  (LIP) 
mitigation  and  blow-in-place  (BIP)  vs.  removal  options.  The  Navy  and  other  DoD  components  have 
responsibility  for  human  safety  and  environmental  stewardship  for  inactive  or  active  ranges  that  have 
been  or  are  currently  being  utilized  for  military  training/readiness  and  for  ordnance  left  underwater  as 
a  result  of  historic  military  activities.  In  addition  to  explosive  blast  (safety)  considerations,  future 
regulatory  emphasis  may  require  an  assessment  of  the  ecological  impact  of  MC  on  the  marine 
environment.  Such  regulations  would  likely  require  an  assessment  of  potential  underwater  ordnance 
contamination  and  mitigation  efforts  that  could  include  water  and  sediment  quality  issues. 

This  effort  directly  addresses  the  amount  of  MC  introduced  into  the  environment  from  the  case  of 
a  single  breached  round  as  a  baseline  scenario.  This  information  can  be  then  used  to  address  the  more 
complex  issues  implicit  for  a  buried  round,  and  also  MC  dispersed  by  a  low-order  detonation. 
Ultimately,  the  information  determined  for  MC  release  from  multiple  breached  shells  can  be  further 
incorporated  into  appropriate  models  predicting  fate  and  transport  of  the  MC  and/or  MC  degradation 
products.  This  study  does  not,  however,  address  the  density,  quantity,  or  distribution  of  munitions  at 
any  particular  site.  Such  data  are  critical  to  assessing  overall  risk  and  must  be  determined  by  other 
methods  or  technologies. 

2.1  SOURCE  CHARACTERIZATION 

Source  characterization  is  the  term  used  to  describe  the  type  and  amount  of  a  contaminant 
introduced  into  the  environment,  generally  as  part  of  an  exposure  pathway  that  requires  the 
assessment  of  risk  to  human  health  and/or  the  environment.  Within  the  context  of  this  proposal  it  is 
defined  as  the  estimation  of  the  mass  of  MC  (TNT,  royal  demolition  explosive;  Cyclotrimethylene- 
trinitramine;  l,3,5-Trinitroperhydro-l,3,5-triazine  (RDX),  high  melting  explosive;  octagen; 
octahydro-l,3,5,7-tetranitro-l,3,5,7-  tetrazocine  (HMX),  and/or  their  breakdown  products  derived 
from  their  military  compositions)  in  both  quantity  and  form  (chemical  species)  that  are  released  from 
the  breached  shell  into  the  water  column  and  sediment. 

To  understand  the  release  of  MC  from  ordnance  items,  it  is  important  to  understand  and  quantify 
the  fundamental  processes  that  govern  the  release  in  the  context  of  the  three  scenarios/conceptual 
models  schematically  shown  in  Figure  1  and  Figure  2. 

These  fundamental  processes  are  described  by  physical  and  chemical  properties  that  ultimately 
define  the  total  magnitude  and  rate  of  MC  release  into  the  surrounding  environment  for  each  of  the 
scenarios  of  interest:  1)  lying  on  top  of  the  sediment,  2)  entirely  buried  in  sediment,  or  3)  as  exposed 
solid  energetic  material  released  directly  to  the  underwater  environment,  e.g.,  in  a  low-order 
detonation.  For  scenarios  1  and  2,  we  assume  that  the  release  through  a  breach  in  a  munition  casing 
can  be  determined  by  the  following  five  key  parameters:  1)  the  start  and  growth  of  the  breach  or  the 
hole  (with  a  size  of  b,  the  radius  of  the  hole);  2)  the  radius  of  the  cavity  formed  due  to  loss  of  mass 
released  from  inside  the  shell,  (R);  3)  the  chemical  property  (dissolution  rate  (p)  from  solid  to 
aqueous  phases  of  the  MC  inside  the  shell  casing);  4)  the  outside  ambient  current  (U)  to  which  the 
casing  hole  is  exposed;  and  5)  mass  of  MC  remaining  inside  (Me).  For  scenario  3,  low-order 
detonation  contamination,  only  parameters  p  (dissolution),  U  (ambient  current),  and  Me  (mass 
remaining)  need  to  be  considered  as  an  extreme  case  where  a  breach  is  infinite  in  size.  It  is  the  goal 
of  this  effort  to  define  and  quantify  the  MC  release  rate  function,  F(b,R,  p,U,Mc),  as  a  function  of  the 
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five  listed  variables,  which  presumably  govern  the  release  rate  of  MC  from  the  munition  casing. 


Figure  1.  Conceptual  model  for  MC  release  from  a  breached  shell  with  a  hole. 

We  have  conducted  both  analytical  and  numerical  modeling  studies  on  MC  release  rates  under 
various  conditions.  To  validate  these  release  rate  results,  an  empirical  study  plan  was  developed  as 
part  of  the  follow-on  (Phase  II  effort  of  this  work).  Results  from  previous  studies  on  the  chemical 
characteristics  of  MCs  in  marine  water  and  sediment  environments  relevant  to  this  study  were 
evaluated  and  used  for  both  the  release  rate  solution  and  subsequent  fate  and  transport  modeling 
study.  For  example,  results  on  chemical  properties,  such  as  solubility  and  dissolution  kinetics  for 
each  of  the  munition  constituents  of  interest,  have  been  used  to  assist  in  defining  the  dependence  of 
the  MC  release  rate  function  (F)  on  the  issues  described  in  Sections  2.1.1  and  2.1.2.. 


Figure  2.  Simplified  MC  release  rate  conceptual  model  configuration. 

2.1.1  Dissolution  Rates  and  Solubility 

Few  measurements  of  aqueous  dissolution  of  MC  from  military  formulations  have  been  conducted. 
The  U.S.  Army  Corps  of  Engineers  (ERDC-EL)  has  shown  (Lynch,  Myers,  Brannon,  and  Delfino, 
2001)  that  solubility  and  dissolution  rates  of  TNT,  RDX,  and  HMX  in  freshwater  were  not 
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significantly  affected  by  pH  over  the  range  4.2  to  6.2,  but  were  significantly  affected  by  temperature 
over  the  range  of  3.1  °C  to  33.3°C.  More  recently,  Lynch  et  al.  (2001)  and  Lynch,  Brannon,  and 
Delfino  (2002a;  2002b)  conducted  studies  measuring  dissolution  rate  with  effects  of  surface  area, 
mixing  rate,  and  temperature  correlated  for  the  separate  dissolution  rates  of  TNT,  RDX,  HMX,  and 
three  military-grade  compositions,  including  Octol  (70/30%  HMX/TNT),  LX-14  (95.5/4.5%  HMX/ 
Estane),  and  Composition  B  (59.5/39.5/1%  RDX/TNT/wax)  (Lynch  et  al,  2001;  Lynch,  Brannon, 
and  Delfino,  2002a;  2002b;  and  Lynch,  Brannon,  Hatfield,  and  Delfino,  2003).  In  those  studies, 
Lynch  et  al  (2002b)  determined  that  dissolution  of  military  compositions  of  explosive  compositions 
occurs  an  order  of  magnitude  slower  than  dissolution  of  the  pure  explosive  compound.  This  is 
probably  due  to  hindered  diffusion  and/or  release  of  explosive  compound  from  the  binding  matrix. 
Increased  salinity  has  been  reported  to  hinder  crystalline  TNT  dissolution,  while  having  little  effect 
on  crystalline  RDX  or  HMX  dissolution,  or  on  dissolution  from  mixtures  within  a  TNT,  RDX,  or 
HMX  military  composition  (Brannon  et  al,  2005).  Environmental  factors  such  as  temperature  and 
pH  can  affect  solubility  and,  therefore,  the  chemical  fate  and  toxicity  (Huang  et  al,  2000). 

Both  dissolution  rates  and  solubilities  of  MCs  in  marine  water  are  used  in  the  release  rate  formulae 
(Brannon  et  al,  2005).  Solid  phase  MC  dissolves  into  the  aqueous  phase  at  a  rate  (dissolution  rate) 
that  competes  with  the  hydrodynamic  diffusivity  created  by  ambient  current  and  the  size  of  the 
breached  hole,  which  act  to  disperse  the  MC  solution  out  of  the  shell. 

2.1.2  Modeling  Fate  and  Transport 

This  study  defines  and  quantifies  the  environmental  process  descriptors  necessary  to  evaluate  the 
fate  and  transport  of  dissolved  MC  in  aquatic  environments.  When  undissolved  MCs  (still  contained 
within  a  shell)  are  introduced  into  the  aquatic  environment,  there  is  no  release  until  the  shell  is 
eventually  breached  by  corrosion  or  other  destructive  processes.  Compounds  inside  are  then  exposed 
to  the  ambient  environment  until  totally  depleted.  Breach  size  {le.,  shell  casing  hole  size)  is  assumed 
to  be  a  function  of  time  in  the  case  of  corrosion.  The  results  of  a  generalized/non-specific  corrosion 
study  initiated  by  Naval  Facilities  Engineering  Service  Center  (Polly,  2012  [in  prep.])  could  be 
integrated  into  the  release  rate  formulae,  in  which  the  hole  size,  b,  will  be  a  time-varying  (growing) 
function  determined  from  the  corrosion  study  result.  Total  MC  mass  remaining  inside  the  shell  is 
another  time-varying  parameter  that  can  be  calculated  by  the  difference  between  the  initial  mass  and 
the  total  MC  mass  released  from  the  shell  from  the  beginning  of  the  formation  of  the  breached  hole. 

Once  released  into  the  environment,  the  MC  sources  undergo  processes  that  govern  fate,  including 
degradation  and  partitioning,  and  transport  processes,  including  diffusivity,  advection,  and  exchanges 
between  the  water  column  and  the  sediment  bed  through  partitioning,  settling,  and  resuspension.  Two 
existing  numerical  models,  the  general  water  quality  model,  WASP  (EPA),  and  the  fine-scaled  fate 
and  transport  model,  TRIM2D  (U.S.  Geological  Survey  and  SSC  Pacific),  are  used  to  simulate  fate 
and  transport  from  the  MC  source  release  from  two  different  locations  in  San  Diego  Bay.  These 
models  have  been  rigorously  scrutinized  and  their  accuracies  demonstrated  by  many  users  nationally 
and  internationally  for  more  than  a  decade.  The  results  of  these  models  should  minimize  uncertainties 
associated  with  model  performance. 

2.2  MODELING  FATE  AND  TRANSPORT  OF  TNT,  RDX,  AND  HMX 

In  Phase  I,  a  conceptual  model  was  first  developed  in  which  a  shell  with  a  circular  breach  hole 
rests  on  the  bottom  with  ambient  current  flowing  tangentially  over  the  hole.  Due  to  the  hydrody¬ 
namic  mixing  through  the  breach  hole,  MC  inside  the  shell  is  dissolved,  an  internal  cavity  is  created 
in  the  solid  MC  matrix,  and  the  MC  solution  is  dispersed  out  of  the  shell  through  the  hole  and  carried 
away  by  the  ambient  current.  Based  on  this  conceptual  model,  the  release  rate  for  a  single  breached 
shell  in  the  marine  environment  was  developed.  Analytically,  the  release  rate  function  can  be 
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expressed  explicitly  with  six  dependent  variables,  which  include  the  ambient  current  speed  (U), 
breach  hole  size  (b),  hydrodynamic  mixing  coefficient  (D),  MC  cavity  radius  inside  the  shell  (R),  the 
MC  dissolution  rate  (p),  and  the  MC  saturation  concentration  (Ds).  Of  these  six  variables,  ambient 
current  and  breach  hole  size  can  be  prescribed  either  from  field  data  or  estimated  by  a  model.  MC 
saturation  concentration  and  dissolution  rate,  both  as  properties  of  the  MC,  can  be  obtained  from 
literature  values  and/or  laboratory  data,  such  as  that  conducted  at  ERDC  (Brannon  et  al,  2005).  The 
hydrodynamic  mixing  coefficient  results  from  interactions  of  the  ambient  flow,  the  internal 
circulation,  and  the  hole  size  and  shape,  and  therefore,  can  be  obtained  by  empirical  data  or  by  a 
numerical  model.  The  semi-analytical  release  rate  function  was  applied,  under  various  hydrodynamic 
and  shell  integrity  conditions,  to  three  common  MCs:  TNT,  RDX,  and  HMX.  For  these  MCs,  the 
time  for  complete  depletion  from  a  breached  shell  under  various  ambient  current  speeds  was 
calculated  to  range  from  less  than  1  year  to  over  1000  years,  depending  on  the  range  of  realistic 
values  selected  for  the  six  variables. 

In  Phase  II,  the  transport  of  MC  release  from  breached  shells  was  studied.  Two  fate  and  transport 
models,  TRIM2D  and  WASP7,  were  evaluated  and  coupled  with  the  MC  release  model  described 
above.  The  WASP7  model  simulates  TNT  concentrations  in  both  the  water  column  and  sediment 
with  box-like  resolutions,  whereas  TRIM2D  provides  simulated  TNT  concentrations  in  the  water 
column  at  much  finer  resolutions.  Both  of  these  models  provide  meaningful  results  that  are  mutually 
consistent  and  complementary,  demonstrating  their  utility  for  developing  an  understanding  of 
munitions  constituent  transport  in  the  marine  environment.  Two  estuarine  scenarios  were  evaluated: 
San  Diego  Bay,  California,  and  the  Elizabeth  River,  Virginia.  Due  to  the  unique  and  distinctive 
characteristics  of  these  two  estuarine  systems,  the  simulated  transport  and  flushing  patterns  of  TNT 
in  these  two  Navy  harbors  exhibit  significance  differences.  Generally,  flushing  in  San  Diego  Bay  is 
driven  by  tides,  whereas  flushing  in  Elizabeth  River  is  driven  primarily  by  the  combined  riverine 
flows  and  tidal  forcing.  It  was  observed  in  the  model  that  differences  in  Total  Suspended  Solids 
(TSS)  settling  velocities  and  associated  TNT  between  the  water  column  and  sediment  bed  is  much 
more  pronounced  for  the  Elizabeth  River  than  for  San  Diego  Bay,  resulting  in  higher  predicted  TNT 
concentrations  in  the  Elizabeth  River  sediments.  Details  of  the  release  rate  formulae  development, 
analysis  and  evaluation  and  partial  verification  using  FLUENT,  empirical  process  descriptors,  and 
preliminary  fate  and  transport  modeling  results  using  TRIM2D  and  WASP7  are  in  Section  4.2  and 
Section  5,  respectively.. 

The  Phase  III  effort  conducted  both  a  numerical  modeling  study  and  an  empirical  laboratory-based 
study  to  verify  the  semi-analytical  release  rate  function.  Results  from  both  the  numerical  and 
empirical  studies  were  compared  with  those  from  the  semi-analytical  release  rate  function. 
Specifically,  the  semi-analytical  release  rate  function  was  compared  and  verified  with  results  from 
the  numerical  computational  fluid  dynamics  (CFD)  model,  FLUENT,  and  also  from  flume-tank 
experiments  conducted  at  the  University  of  Wisconsin-Milwaukee.  This  report  documents  study 
results  for  1)  effect  of  opening  shapes,  2)  internal  circulation,  3)  probabilistic  modeling  of  the  release 
rate,  and  4)  an  empirical  study  for  unburied  shells  (Figure  3). 

Numerical  simulations  were  conducted  to  1)  verify  the  release  rate  function,  and  2)  investigate  the 
effect  on  the  release  rate  from  different  shapes  of  the  hole  opening.  Mass  release  rates  were 
calculated  for  breach  holes  with  different  shapes,  with  the  same  opening  areas,  including  a  circle,  a 
triangle,  and  two  rectangles  with  aspect  ratios  of  4:1  and  10:1,  respectively.  The  high  aspect  ratio 
rectangular  shape  was  investigated  for  simulating  cracks  in  a  shell. 

The  variables  for  the  mass  release  rate  function  can  be  treated  in  both  a  deterministic  and 
probabilistic  fashion.  The  probabilistic  modeling  has  two  components.  The  first  addresses  real-world 
scenarios  and  focuses  on  the  variability  associated  with  selected  input  parameters,  and  the  second  is 
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focused  on  (defined  or  undefined)  uncertainties  associated  with  those  input  parameters.  For  the  first 
scenario,  ambient  current  amplitudes  and  hole  sizes  were  assigned  probabilistic  distribution 
functions,  then  Monte-Carlo  simulations  were  conducted.  Predicted  mass  release  rates  for  these 
scenarios  were  associated  with  probabilistic  distributions.  Subsequently,  the  Monte-Carlo-based 
probabilistic  model  allowed  for  uncertainty  and  sensitivity  analyses  of  the  predicted  mass  release 
rates. 

The  empirical  study  was  conducted  at  the  flume  tank  of  the  University  of  Wisconsin- 
Milwaukee.  The  tank  was  designed  to  run  flows  with  a  current  speed  range  of  ~  7  to  45  em/s. 
Semi-cireular  hollow  domes  (purchased)  with  different  radii  were  sequentially  glued  to  an 
acrylie  sheet.  Holes  of  different  radii  were  drilled  at  the  center  of  an  acrylie  sheet  forming  the 
impermeable  closure  to  the  top  of  the  half  dome.  Fluoreseein  was  used  as  the  surrogate  dye  for 
MC  and  was  plaeed  inside  the  half  dome,  which  was  placed  under  the  acrylie  sheet  over  which 
tank  water  flowed  at  the  designed  speed.  Particle  image  veloeimetry  (PIV)  was  used  to  measure 
flow  speed  inside  and  outside  the  half  dome.  Released  dye  from  the  hole  was  measured  ~20  em 
downstream  by  the  laser-induced  fluorometry  (LIF)  device.  Both  data  were  measured  as  images, 
and  image  proeessing  algorithm  software  was  used  to  ealculate  flow  speed  from  PIV  data  and 
dye  eoneentrations  from  the  LIF  data. 


Figure  3.  Study  plan  and  tasks. 

2.3  MC  RELEASE  AND  TRANSPORT  BURIED  IN  SAND 

Pore  water  flows  in  coastal  water  are  affected  by  multiple  factors,  including  advection  induced  by 
surface  waves,  ground  water  seepage,  and  mixing  by  both  flow-induced  mixing  and  molecular 
diffusivity.  In  near-shore,  coastal  regions,  surface  waves  exert  pressure  on  the  sediment  bed,  and 
hence,  pore  water  in  the  sediment  bed  is  subject  to  wave  actions.  Such  actions  and  their  impacts  on 
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pore-water  advection,  although  minute,  are  continuous,  and  thus  the  accumulated  effects  on  the 
release  rate  are  likely  not  negligible. 

Similar  to  wave  actions,  groundwater  flows  in  coastal  and  near-shore  regions  have  effects  on  the 
advection  and  mixing  of  MC  release.  Groundwater  seepage  in  these  regions  results  from  multiple 
sources.  For  all  of  these  forcing  functions,  the  seepage  flows  follow  Darcy’s  law,  which  correlates 
flow  with  pressure  (hydraulic  head)  gradients. 

2.3.1  Groundwater  Seepage 

Groundwater  in  coastal  region  sediments  originates  from  the  interactions  of  ocean  water  and  inland 
surface  waters  (rivers,  streams,  lakes)  and  groundwater.  Generally,  ocean  water  and  inland  water 
meet  and  interact  near  the  shoreline  and  beach  areas.  These  interfaces  are  composed  of  porous  media 
through  which  interactions  and  exchanges  of  inland  freshwater  with  the  saline  ocean  water  take 
place.  For  a  cross-section  view  cubing  across  the  shoreline  (Figure  4),  the  direction  of  pore-water 
flow  is  seaward  under  low  tide  and  high  river  flow  or  high  water  table  conditions.  Conversely,  pore- 
water  flow  is  landward  under  high  tide,  low  river  flow,  or  low  water  table  conditions.  Pore-water 
advection  in  these  coastal  and  shoreline  regions  has  been  studied  extensively  in  both  field  efforts  and 
modeling  studies.  In  general,  a  groundwater  model,  such  as  MODFLOW  or  SUTRA,  can  be  applied 
to  quantify  and  predict  dynamics  of  groundwater  and  pore-water  flows.  Based  on  the  fundamental 
underlying  principles  of  these  models,  a  simple  two-dimensional  (2-D)  model  has  been  developed  for 
predicting  the  pore-water  advection  in  the  sediment  flume  laboratory  study. 


Figure  4.  Interactions  between  groundwater  and  surface  ocean  water. 

2.3.2  Pore-water  Advection  Induced  by  Wave  Effect 

Pore- water  velocities  associated  with  current-induced  pressure  gradients  (Figure  5)  are  typically  in 
the  range  of  several  centimeters  (0.1  to  9.2  cm)  per  hour  (Huettel  and  Gust,  1992;  Huettel,  Forster, 
Kloser,  and  Fossing,  1996).  Advective  pore-water  flow  in  the  natural  environment  has  also  been 
reported.  Tracing  of  dye  plumes  injected  into  the  sand  at  shallow,  rippled  Mediterranean  sites 
revealed  the  presence  of  pore-water  advection  in  the  presence  of  wave-induced  oscillating  bohom 
flow  (Precht  and  Huehel,  2004).  In  Mid-Atlantic  Bight  sediments  that  were  exposed  to  combined 
tidal  and  wave-driven  currents,  the  presence  of  fast  vertical  pore-water  transport  became  apparent 
from  the  short  time  it  took  iodide,  released  above  the  sediment  surface,  to  arrive  at  the  tip  of  a 
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voltammetric  sensor  20  mm  below  the  sediment  surface  (Reimers  et  al,  2004).  Even  relatively  slow 
advective  pore- water  circulation  at  velocities  of  a  few  centimeters  per  hour  resulted  in  intense 
interfacial  exchange  of  solutes.  In  flume  experiments  with  stained  pore  waters,  Huettel  and  Gust 
(1992)  observed  interfacial  pore-water  exchange  in  the  presence  of  small  mounds  (height  <  25  mm) 
and  quantified  the  rates  in  terms  of  dye  release  to  the  overlying  water. 

In  a  low-pressure,  bottom-sea  mound  (sand  ripple,  as  per  Figure  5),  pore  water  was  transported 
upwards  from  >100  mm  sediment  depth  and  released  at  the  sediment  surface.  Near  the  mounds,  dye 
release  exceeded  diffusion  by  three  orders  of  magnitude.  Based  on  these  measurements,  Huettel  and 
Gust  (1992)  estimated  that  advection  induced  by  mounds  at  natural  densities  in  the  North  Sea 
intertidal  would  affect  the  entire  surface  area  and  result  in  a  fold  increase  of  fluxes  when  compared  to 
diffusion.  Using  a  similar  experimental  approach,  other  researchers  (Precht  and  Huettel,  2003) 
quantified  solute  exchange  across  rippled  beds  in  a  wave  tank  under  oscillating  currents  and  found 
release  rates  that  were  50  times  greater  than  diffusive  exchange. 

2.3.3  Wave  Actions 

Coastal  waves  impose  pressure  on  the  sediment  bed  in  coastal  shallow  regions.  Due  to  their 
dispersive  nature,  ocean  waves  propagate  at  various  speeds  according  to  their  wavelengths  and 
frequencies.  When  waves  propagate  from  deep  ocean  water  toward  the  shallow  coastal  water,  they 
gradually  impose  increasing  pressure  onto  the  sediment  bed.  As  these  waves  are  restricted  by  the 
shallow  near-shore  regions,  the  wave  heights  grow  and  eventually  break  in  the  surf  zone,  which 
further  impose  impacts  to  the  bottom  sediments.  These  wave-induced  impacts  on  sediment  and  pore 
water  are  briefly  discussed  here. 

Generally,  ocean  waves  can  be  approximated  by  the  sum  of  multiple  trains  of  single  linear  waves, 
each  of  which  has  an  amplitude  expressed  as 

T]{x,t)  =  Hsm{kx-cot) ,  (1) 

where  H  is  a  constant  representing  the  maximum  amplitude  of  the  wave  train.  The  wave  number, 
k  =  InlX,  and  the  frequency,  oo,  are  related  by  the  dispersion  equation: 

0)^  =  gk  tanh(  kh) 

=  gk  for  deep  water;  h/A  >  1/4 
=  gk^h  for  shallow  water;  h/A  <1/11, 


where  g  is  the  gravity  constant  (9.8  m/s^),  X  is  the  wave  length,  h  is  the  local  depth,  and  tanh  is  the 
hyperbolic  tangent  function. 

The  speed  of  the  traveling  wave  (traveling  velocity)  can  be  expressed  as  the  ratio  of  the 
wavelength  over  the  period: 
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where  T  is  the  wave  period  (in  seconds). 

As  can  be  shown  in  Eqs.  (2)  and  (3),  deep-water  waves  are  dispersive  long  waves,  with  low 
frequencies  traveling  faster  than  short  waves  with  high  frequencies.  Note  that  in  the  equations  both 
deep-  and  shallow-water  waves  are  defined  by  the  parameter,  kh,  which  is  a  ratio  of  water  depth 
versus  wave  length.  Therefore,  a  wave  with  the  same  wavelength  is  defined  as  a  shallow  water  or 
deep  water  wave  according  to  the  parameter,  kh.  In  addition,  for  a  group  of  waves  originated 
offshore,  long  waves  travel  faster  and  reach  the  surf  zone  first.  As  these  long  deep-water  waves  travel 
to  the  shallower  regions,  they  gradually  become  shallow-water  waves,  with  traveling  velocity 
restricted  by  the  local  depth,  again  as  shown  in  Eq.  (3). 

Pressure  exerted  by  a  wave  train  to  the  bottom  sediment  can  be  expressed  as 
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The  pressure  exerted  by  the  wave  train  on  the  bottom  propagates  along  the  x-direction  and  is 
oscillatory  in  time.  Pressure,  and  its  gradient  exerted  by  the  wave  at  the  sediment  water  interface, 
also  decays  exponentially  with  depth.  Therefore,  surface  waves  lose  their  effects  rapidly  as  water 
depth  increases.  Eq.  (4)  provides  the  boundary  condition  that  generates  pore-water  flows,  which  is 
also  oscillatory.  In  reality,  advection  in  pore  water,  resulting  from  the  pressure  differences  of  the 
waves,  dissipates  fast  inside  the  sediment  bed,  because  ocean  waves  are  random  with  multiple  wave 
directions,  which  counteract  wave  actions  from  multiple  directions  and  wave  components. 

2.3.4  MC  Release  from  the  Deep  Sea  Environment 

As  discussed  above,  groundwater  seepage  and  wave  actions  can  produce  pore-water  advection  in 
shallow,  coastal  regions,  where  the  interactions  of  fresh  surface  water  and  saline  ocean  water  take 
place  and  surface  wave  actions  can  reach  the  bottom.  These  effects  are  reduced  significantly  (expo¬ 
nentially)  in  deep-sea  environments.  However,  MC  release  rate  is  affected  by  deep-sea  currents, 
which  are  affected  by  global  circulation,  and  altered  dissolution  rate  due  to  low  water  temperature. 

In  addition  to  very  low  water  temperature,  the  deep  sea  is  also  characterized  by  high  dissolved 
oxygen,  darkness  with  no  sunlight,  and  currents  with  variable  magnitudes.  All  of  these  features  have 
potential  impacts  on  the  release  rates  of  MC  in  deep-ocean  regions.  To  study  MC  release  from 
breached  shells  in  deep-sea  environment,  the  semi-analytical  release  function  developed  and 
validated  in  this  study  is  further  analyzed  by  evaluating  the  relative  magnitudes  for  three  of  the  five 
governing  variables,  including  1)  ambient  current  (U),  2)  hydrodynamic  mixing  coefficient  (D),  and 
3)  dissolution  rate  constant  (p).  The  other  two  variables,  hole  size  and  internal  cavity  volume,  would 
have  no  significant  differences  between  deep  and  shallow  seawater  scenarios. 
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Deep-Sea  Currents 

Deep-sea  currents  are  primarily  driven  by  four  major  forcing  components,  including  energy  from 
the  sun,  the  rotation  of  the  earth,  the  shape  of  land  masses,  and  shape  of  the  ocean  floor.  Solar  energy 
is  more  intense  at  the  equator  than  at  the  poles,  causing  the  air  and  water  near  the  equator  to  be  hotter 
than  at  higher  latitudes.  The  heated  warm  water  near  the  equator  moves  the  water  along  the  density 
gradient  northward  (Figure  6).  For  example,  the  Gulf  Stream  is  driven  by  the  warm  water  from  the 
equator  traveling  north  along  the  east  coast  of  America  (Gaskell,  1973).  As  the  current  travels  further 
north,  one  section  moves  south  down  the  west  coast  of  Africa,  before  completing  a  full  circle  to 
return  to  the  Gulf  of  Mexico.  The  second,  often  referred  to  as  the  North  Atlantic  Drift,  moves  further 
north  up  the  west  coast  of  Ireland  and  Scandinavia  and  cools  at  high  altitudes  and  then  sinks  as  a 
deep-water  cold  current  back  along  the  east  coast  of  America.  Due  to  the  rotation  of  the  earth,  the 
Coriolis  Effect  and  trade  winds  force  the  currents  to  flow  in  a  clockwise  direction  (Gross,  1982).  In 
addition  to  the  above,  land  masses  and  their  distributions  also  dictate  the  currents  and  circulations  in 
waters  bounded  by  the  land. 

Magnitudes  of  deep-sea  currents  have  been  reported  to  vary  with  geographical  location,  ranging 
from  2  to  10  cm/sec  on  the  low  end  (Adams,  1999a  and  1999b),  to  4  to  14  cm/sec  (Gross,  1990),  and 
to  the  high  end  of  -150  cm/sec  (Coble,  Murray,  and  Rice,  1987),  within  the  Gulf  Stream.  These 
current  speed  ranges  are  equivalent  to  those  for  tidal  currents  in  estuaries,  coastal  regions,  and 
harbors  for  which  the  MC  release  study  was  designed.  Therefore,  impacts  to  the  MC  release  rates 
from  deep-sea  currents  depend,  to  a  great  extent,  on  the  magnitude  of  the  currents,  in  a  pattern 
similar  to  those  in  the  shallow  estuaries,  harbors,  and  coastal  regions. 


Low  Pressure 


Figure  5.  Schematic  of  pore-water  flow  field  driven  by  oscillating  flow  interacting  with  a  sediment 
ripple.  Top  and  center  show  the  pore-water  flow  field  that  would  develop  under  steady  unidirectional 
flow  in  opposing  directions;  the  bottom  drawing  shows  the  averaged  pore-water  flow  field  under 
oscillating  flow  as  induced  by  surface  gravity  waves  (from  Precht  and  Huettel,  2003). 
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Hydrodynamic  Mixing 

Hydrodynamic  mixing  also  plays  a  role  in  mixing  and  exchanging  solute  with  ambient  water 
through  the  hole.  It  is  also  demonstrated  that  the  hydrodynamic  mixing  is  strongly  correlated  with  the 
magnitude  of  the  ambient  current  and  size  of  the  internal  cavity,  and  to  a  lesser  degree,  with  the 
breach  hole  size.  Since  the  internal  cavity  size  and  its  growth  are  directly  related  to  the  overall 
release  rate,  i.e.,  the  reduction  in  munition  fill  through  MC  dissolution  and  release,  the  internal  cavity 
factor  can  be  neglected  since  it  is  a  result  of  the  overall  release  rate.  In  the  deep  sea,  wave  effects  and 
groundwater  are  negligible,  which  leaves  deep-sea  current  as  the  single  factor  that  determines  the 
hydrodynamic  mixing  parameter,  which  follows  the  release  rates  validated  in  this  study. 


Figure  6.  Global  surface  current  system  under  average  conditions  for  winter  months  in  the  Northern 
Hemisphere.  Warm  currents  are  shown  as  solid  red  arrows,  and  cold  currents  as  dashed  blue 
arrows. 

Hydrodynamic  Mixing 

Hydrodynamic  mixing  also  plays  a  role  in  mixing  and  exchanging  solute  with  ambient  water 
through  the  hole.  It  is  also  demonstrated  that  the  hydrodynamic  mixing  is  strongly  correlated  with  the 
magnitude  of  the  ambient  current  and  size  of  the  internal  cavity,  and  to  a  lesser  degree,  with  the 
breach  hole  size.  Since  the  internal  cavity  size  and  its  growth  are  directly  related  to  the  overall 
release  rate,  i.e.,  the  reduction  in  munition  fill  through  MC  dissolution  and  release,  the  internal  cavity 
factor  can  be  neglected  since  it  is  a  result  of  the  overall  release  rate.  In  the  deep  sea,  wave  effects  and 
groundwater  are  negligible,  which  leaves  deep-sea  current  as  the  single  factor  that  determines  the 
hydrodynamic  mixing  parameter,  which  follows  the  release  rates  validated  in  this  study. 

Dissolution  Rate  Constants  (^) 

Dissolution  rate  constants  are  temperature  dependent,  which  may  be  most  influenced  by  the  low 
water  temperatures  found  in  deep-sea  environments.  Before  we  investigate  the  temperature  effect  on 
the  release  rates,  temperature  in  the  deep-sea  environment  is  briefly  discussed. 

Temperature  of  ocean  water  is  influenced  by  multiple  factors,  including  heat  from  the  sun,  salinity, 
and  water  mass  transport  (circulation  and  the  mixing).  Much  of  the  ocean  water  (over  90%)  is  found 
below  the  thermocline,  which  is  located  50  to  100  m  in  the  deep  ocean.  Temperature  should  have  an 
effect  on  the  dissolution  rate,  which  is  a  parameter  of  the  release  rate  function.  Dissolution  rate 
constants  for  TNT,  RDX,  and  HMX  in  freshwater  have  been  published  for  freshwater  (Lynch, 
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Brannon,  and  Delfino,  2002a)  and  for  seawater  (Brannon,  et  al.,  2005).  Dissolution  rates  for  the  three 
munition  constituents,  including  TNT,  HMX,  and  RDX,  do  not  vary  by  much  when  compared  to 
freshwater  and  seawater  (Figure  7,  Brannon,  et  al.,  2005)  The  rate  constants  for  freshwater  (Lynch, 
Brannon,  and  Delfino,  2002a)  were  measured  at  three  different  temperatures,  10°C,  20°C  and  30°C, 
respectively.  Generally,  temperatures  for  typical  coastal  surface  water  are  within  these  ranges  and 
measured  rate  constants  can  be  directly  used  for  these  regions.  However,  in  the  deep  ocean,  tempera¬ 
ture  can  drop  rapidly  below  the  thermocline  to  the  range  of  0  to  4°C.  Dissolution  rate  constants 
measured  for  the  normal  shallow  water  temperatures,  namely  10°C,  20°C,  and  30°C,  can  be  linearly 
extrapolated  for  low  temperature  conditions. 


Figure  7.  Dissolution  rates  of  TNT,  RDX,  and  HMX  from  TNT,  military  formulations  of  Composition  B, 
and  Octol  (Brannon  et  al.,  2005). 

Figure  8  shows  the  dissolution  rates  of  TNT,  RDX,  and  HMX,  extrapolated  for  4°C  water  from 
those  measured  by  Lynch,  Brannon,  and  Delfino  (2002a).  Compared  to  the  20°C  coastal  surface 
water,  the  dissolution  rates  are  reduced  by  a  factor  of  2.7  for  TNT,  a  factor  of  10  for  RDX,  and  a 
factor  of  7  for  HMX.  The  temperature  effect  is  much  greater  than  the  salinity  effect.  This  is 
significant  because  the  rate  constants  directly  govern  the  dissolution  rate  in  a  linear  fashion  and  the 
reduced  dissolution  in  a  cold,  deep-water  environment  dictates  that  it  should  take  a  longer  time  for 
MC  depletion  from  breached  shells  in  these  environments. 
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Figure  8.  Dissolution  rates  for  freshwater  as  a  function  of  temperature  (symbols  with  solid  lines  are 
empirical  freshwater  data  from  Lynch,  Brannon,  and  Delfino  (2002a).  Dashed  line  represents  best 
fits  of  the  data  for  extrapolation  to  low  temperature. 
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3.  MATERIALS  AND  METHODS 


3.1  RELEASE  FUNCTION  OF  MC  FROM  BREACHED  SHELLS 

The  approach  to  estimate  the  release  rate  of  unexploded  munition  constituents  (MCs)  from 
breached  shells  into  estuarine  and  coastal  waters  has  included  the  development  of  an  analytical 
release  function  and  partial  verification  using  a  numerical  model,  FLUENT.  The  analytical  release 
function  is  derived  with  the  assumption  that  MC  is  released  through  the  breached  hole  and  exposed 
to  ambient  hydrodynamic  conditions.  Specifically,  the  MC  release  through  a  breach  in  a  munition 
casing  is  determined  by  the  following  five  key  parameters:  1)  the  start  and  growth  of  the  breach  or 
the  hole  in  the  shell  casing  (with  a  size  of  b,  the  radius  of  the  hole),  2)  the  radius  of  the  cavity  formed 
due  to  loss  of  mass  released  from  inside  the  shell  (R ),  3)  the  chemical  property  (dissolution  rate  (p) 
from  solid  to  aqueous  phases  of  the  MC  inside  the  shell  casing),  4)  the  ambient  current  (U)  to  which 
the  casing  hole  is  exposed,  and  5)  mass  of  MC  remaining  inside  (Me).  This  effort’s  goal  is  to  define 
and  quantify  the  release  function,  F(b,  R,  p,  U,  Me),  as  a  function  of  the  five  variables,  which 
presumably  govern  the  release  rate  of  MC  from  the  munition  casing 

This  analytical  release  function  was  validated  by  the  numerical  model,  FLUENT,  which  is  a 
commercial  computational  fluid  dynamic  (CFD)  software  program.  The  munition  source  is 
considered  a  passive  scalar  concentration  governed  by  a  fluid  transport  equation,  which  is  solved 
concurrently  with  flow  equations.  The  MC  release  rate  is  calculated  from  simulation  results.  The 
dependency  of  the  release  rate  on  the  breach  size,  ambient  current,  and  shell  thickness  were  exam¬ 
ined  with  FLUENT  modeling.  Based  on  simulation  results,  this  analytical  equation  for  the  release 
rate  is  proposed,  which  can  be  easily  implemented  as  a  source  term  into  any  fate  and  transport  model. 

3.2  DISPERSION  OF  MC  FROM  BREACHED  SHELLS 

In  the  MC  release  problem,  it  is  assumed  that  MC  is  released  from  the  shells  through  a  breach  or 
hole  in  the  shell  casing,  which  is  exposed  to  ambient  currents.  The  size  of  the  hole  is  assumed  to  be 
small  compared  to  the  overall  size  of  most  shells,  with  curvatures  that  are  small  enough  not  to  affect 
the  ambient  currents  at  the  outer  surface  of  the  shell  as  in  the  conceptual  model  (Figure  1).  Under 
these  assumptions,  the  MC  release  process  can  be  approximated  by  a  flat  plate  as  in  the  simplified 
conceptual  model.  Figure  2,  in  which  MC  is  released  from  a  hole  on  the  plate  over  which  ambient 
current  flows. 

It  is  assumed  that  most  of  the  total  release  occurs  during  a  stage  when  a  relatively  large  cavity 
(compared  to  the  breach  hole  or  crack)  has  formed  due  to  dissolving  of  the  MC  matrix.  Thus,  the 
overall  release  rate  is  governed  by  the  advection  and  diffusion  of  MC  concentration  within  the  inner 
cavity,  which  is  in  turn  determined  by  the  dissolution  of  MC  at  the  solid-solution  interface  and  the 
exchange  flow  at  the  breach  (boundary  conditions).  The  driving  force  for  MC  leaching  from  the  inner 
cavity  is  the  hydrodynamic  diffusion  caused  by  ambient  current  inside  the  cavity,  which  generates 
inner  circulation.  We  also  assume  that  the  local  shear  stress  near  the  breach  hole  largely  determines 
the  inner  circulation  of  the  MC  solution,  e.g.,  different  external  flow  conditions  with  similar  near¬ 
breach/hole  hydrodynamic  conditions  result  in  a  similar  inner  circulation  pattern,  resulting  in 
approximately  the  same  mass  depletion  rate. 

Based  on  the  above  assumptions,  the  MC  flux  can  be  examined  with  a  simplified  conceptual  model 
consisting  of  a  uniform  ambient  fluid  flowing  over  a  flat  plate  that  represents  the  surface  of  a  shell  as 
in  Figure  1  and  Figure  2.  The  shell  breach  is  modeled  as  a  circular  hole  with  radius  h.  Water  flows 
over  the  hole  with  speed  U.  The  typical  flow  inside  the  cavity  is  characterized  with  three-dimensional 
(3-D)  structures.  However,  we  assume  a  cavity  with  a  spherical  shape  with  radius  R,  i.e.,  diffusion/ 
mixing  process  (by  molecular  diffusion  and  flow-induced  stirring)  is  considered  homogeneous.  The 
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3-D  diffusion/mixing  problem  inside  the  casing  is  thus  simplified  to  a  one-dimensional  (1-D) 
process.  The  controlling  equation  can  be  written  as 


dC  d  (  2dC'] 

—  =  D— ^ —  r  — 
6t  5r  dr ) 


(5) 


where  C  is  concentration  as  a  function  of  radial  direction,  r.  and  time,  D,  is  a  dispersion/diffusivi-  ty 
coefficient,  which  may  include  the  combined  effects  of  molecular  diffusion.  Dm,  and  circulation 
induced  mixing.  Da,  {i.e.,  D  =  Dm+Da),  where  Dm  is  known  as  a  physical  property  of  seawater  and 
MC  and  Da  is  an  hydrodynamic  function  induced  by  ambient  current  speed,  U,  which  will  be 
obtained  empirically  from  the  results  of  the  numerical  model,  FLUENT. 


At  an  equilibrium  state,  the  general  solution  to  Eq.  (5)  is 


C{r)  =  B - , 

r 


where  r  is  hole  radius,  and  A  and  B  are  constants. 


The  boundary  condition  at  the  solid-liquid  interface  is  prescribed  as 


(6) 


-D— =  //(C,-C)  @r  =  R,  (7) 

dr 

where  //  is  the  dissolution  speed  [cm/s],  which  is  obtained  by  dividing  the  dissolution  rate  [mg/cm^/s] 
over  the  saturation  concentration  Cs  (solubility),  [mg/L].  This  boundary  condition  reflects  the  fact 
that  flux  by  diffusion  at  the  solid/liquid  interface  is  in  equilibrium  with  the  dissolution  process. 

At  the  breach  hole,  the  flux  of  C  caused  by  mean  current  is  considered  to  be  linearly  proportional 
to  U,  the  local  C  and  the  2-D  area  of  the  hole,  i.e.,  MC  flux  is  prescribed  as  aUCnb^,  which  is  in 
equilibrium  with  diffusive  flux  from  just  inside  the  hole.  The  coefficient  a  (approximately  equal  to  1) 
is  a  model  parameter  that  is  a  function  of  the  local  geometry  of  the  breach  hole/crack  and  the  thick¬ 
ness  of  the  shell.  Thus,  we  can  write  the  boundary  condition  at  the  breach  hole  as 

Ijih^  D - =  aUC7ib^  @r  =  b 

dr 


or 

dC  _aU  ^ 
dr  2D 


@r  =  b . 


(8) 


Substituting  boundary  conditions,  Eqs.  (7)  and  (8),  into  the  general  solution,  Eq.  (6),  we  have 
obtained  the  following  concentration  field  inside  the  cavity; 


C(r)  = 


_ _ f  2D  J._n 

2.D  ^  _i_  1  1  y^aUb^  b  r) 

aUb^  b  R 


(9) 
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The  MC  release  rate  can  be  obtained,  by  definition,  as  follows; 


F  =  aUC{b)7ib^ 


2D  Dll' 

aUb^  juR^  ly  R 


(10a) 


F  is  the  release  rate  function,  which,  as  depicted  in  Eq.  (10a),  is  a  closed- form  solution  with  the 
five  variables,  including  hydrodynamic  diffusivity  coefficient  (D),  current  (U),  hole  size  (b),  cavity 
radius  (R)  and  dissolution  speed  of  MC  from  solid  to  aqueous  phase  (p).  The  model  parameter  a  is 
defined  as  a  geometry  factor  (Eq.  (8)). 

Eq.  (10a)  can  also  be  expressed  in  a  non-dimensional  form  as  follows: 
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The  denominator  of  Eq  (10b)  includes  three  non-dimensional  parameters,  including  the  current- 
based  Reynolds  Number  (Ub/D);  the  model  parameter  a  and  the  number  2,  which  are  both  constants; 
the  dissolution-based  Reynolds  Number  (pR/D);  and  the  hole-to-cavity  size  ratio  (b/R).  The 
hydrodynamic  properties,  D  and  U,  can  be  described  by  field  data  or  a  numerical  fluid  dynamic 
model.  The  coefficient  k  in  the  numerator  can  be  obtained  empirically  either  by  laboratory 
experiment  or  from  the  FLUENT  model.  In  our  numerical  modeling  experiments  using  FLUENT,  to 
be  discussed  later,  the  value  of  k  is  between  0.015  and  0.035,  and  k  =  0.023  if  shell  thickness  is 
neglected. 

Generally,  the  breached  hole  is  much  smaller  than  the  cavity  during  scenarios  where  mass  flux  is 
most  significant.  As  a  result,  the  1/R  term  in  the  denominator  of  Eq.  (10a)  can  be  neglected,  which 
will  help  to  highlight  the  behavior  of  the  release  rate  function  for  the  two  extreme  conditions 
described  in  Sections  3.2.1  and  3.2.2. 

3.2.1  Current-Controlled  Release  Function  (pR^  »  aUb^) 

When  the  dissolution  speed  (p)  is  much  larger  than  the  ambient  current  (U)  or  the  cavity  (R)  is 
much  larger  than  the  hole  (b),  MC  solutions  inside  the  cavity  are  near  saturation  and  the  release  rate 
of  MC  is  controlled  by  the  ambient  current  speed  U  and  the  hole  size.  The  release  rate  function  can 
be  simplified  explicitly  as 


_  27:DC^ 

2D  ^  r  (11) 

aUb^  b 
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The  MC  release  rate  increases  with  ambient  current  speed,  the  area  of  the  hole,  and  the  saturation 
concentration  of  the  MC. 

3.2.2  Dissolution-Rate-Controlled  Release  Function  (pR^  «  aUb^) 

For  this  case,  release  of  MC  is  controlled  by  the  dissolution  speed  (rate),  since  ambient  current 
and/or  the  hole  is  so  large  that  MC  solution  dissolved  from  the  solid  phase  is  dispersed  out  of  the 
breached  shell  immediately; 


D  V 

(12) 

« ;  if  «  — 

^  D  R 

When  the  hole  size  is  relatively  small  (compared  to  the  parameter,  pR^/D),  the  hole  has  a  damping 
effect,  restricting  the  release  rate  of  MC  out  of  the  shell.  When  the  hole  size  becomes  relatively  large 
(compared  to  pR^/D),  the  damping  effect  of  the  hole  diminishes  and  the  release  of  MC  is  directly 
proportional  to  the  surface  area  of  the  cavity  (27iR^)  and  the  mass  dissolution  rate  (pCs),  as  indicated 
inEq.  (12). 

3.2.3  Low-Order  Detonation  Scenario  Release 

In  the  case  where  shells  are  broken  apart  and  MC  is  fully  exposed  to  the  ambient  current,  MC 
dissolution  from  the  solid  phase  to  the  aqueous  phase  is  controlled  by  the  dissolution  rate.  Once 
dissolved,  the  MC  plume  is  advected  away  by  the  ambient  current.  Therefore,  under  this  low-order 
detonation  scenario,  MC  release  follows  the  dissolution-controlled  release  process.  The  release 
function,  Eq.  (12),  derived  for  release  through  a  hole,  can  be  modified  for  this  scenario: 

F  =  As/jCs,  (13) 

where  the  term  1/b  in  Eq.  (12)  is  infinitely  small  (hole  is  so  large  that  1/b  ~  0)  and  the  solid- water 
interfacial  area  of  the  semi-sphere,  27iR^,  is  replaced  by  surface  area  of  the  MC  solid-water  interface 
(As).  For  this  scenario,  ^^^may  be  obtained  from  field  observation/data.  Therefore,  release  of  MC 
under  a  low-order  detonation  scenario  depends  on  the  integrity  of  the  exposed  MC,  represented  by 
the  total  contact  surface  area  between  the  remaining  active  MC  matrix  and  the  aqueous  phase. 

3.3  EMPIRICAL  STUDY  -  UNBURIED  SHELL 
3.3.1  Experiment  Test  Platform 

All  mass  release  experiments  were  conducted  in  a  tiltable  re-circulating  open-channel  flume  in  the 
hydraulic  lab  at  University  of  Wisconsin-Milwaukee.  The  flume  consists  of  a  head  tank,  a  tail  tank, 
and  a  test  section  (Figure  9  and  Figure  10).  The  channel  is  60-feet  long,  and  2.5-feet  wide,  with  a 
water  depth  up  to  2  feet.  The  two  side  walls  of  the  test  section  are  made  of  glass,  allowing  excellent 
optical  viewing.  The  slope  of  the  channel  is  adjustable,  with  the  water  flow  driven  by  a  centrifugal 
pump.  The  flow  rate  of  the  channel  can  be  adjusted  by  varying  the  pump  speed,  the  channel  slope, 
and  using  a  valve  simultaneously.  Water  is  sucked  from  the  tail  tank  by  the  pump,  carried  through  a 
pipe  underneath  the  flume,  and  upwells  into  the  head  tank.  It  then  flows  through  a  series  of  grids  that 
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regulate  the  flow  condition  to  generate  a  uniform  mean  current  with  a  homogeneous  turbulence 
distribution  in  the  test  section.  All  measurements  were  conducted  at  a  location  about  35  feet 
downstream  of  the  inlet,  where  a  mass  release  test  platform  was  placed. 


Flush  mounted 


Figure  9.  Re-circulating  flume  tank  and  the  conceptual  release  from  the  half  dome. 


Figure  10.  Test  section  of  open  channel  flume  at  UWM  hydraulic  lab. 
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A  mass  release  platform  was  constructed  to  create  an  ideal  uniform  flow  over  a  smooth  surface. 
The  platform  was  made  of  polycarbonate  plates  with  thickness  =  1  inch,  a  total  length  of  4  m,  and 
width  of  0.75  m  (about  equal  to  the  width  of  the  flume).  The  platform  was  placed  in  the  open 
channel,  supported  by  aluminum  rods  and  was  elevated  to  10  cm  above  the  flume  bed.  Breached 
cavities  were  modeled  from  semispherical  acrylic  domes.  Breached  surfaces  were  made  of  1/8-inch 
thick  circular  acrylic  sheets,  which  were  glued  to  the  cavity  dome  using  silicone  sealer.  Circular 
holes  or  line  cracks  of  various  sized  were  machined  by  drilling  through  the  covering  sheets  over  the 
glued  cavity  domes.  Modeled  cavities  with  breached  holes/cracks  were  then  flush-mounted  onto  the 
flow  platform  with  the  center  of  holes/cracks  located  3  meters  downstream  of  the  leading  edge  of  the 
platform.  Figure  1 1  shows  the  designed  interchangeable  mass  release  test  platform.  The  experimental 
system  was  designed  in  such  a  way  that  cavities  with  different  breach  sizes,  shapes,  and  inner  radii 
could  be  easily  mounted  on  the  same  platform.  For  each  designed  flow  condition,  experiments 
commenced  once  the  flow  reached  a  uniform  condition  along  the  mean  stream  direction  at  the 
location  of  the  breach. 


Figure  1 1 .  Model  flow  platform  and  shell  cavity. 

3.3.2  Measurement  Techniques 

To  characterize  the  external  shearing  flow  above  the  breach,  as  well  as  the  induced  in-cavity 
circulation,  PIV  was  applied  to  measure  the  2-D  instantaneous  velocity  field.  The  PIV  system 
consisted  of  a  high-resolution,  cross-correlation  CCD  camera,  a  dual-head  Q-switched  laser,  a  set  of 
optics,  and  a  PC -based  image  acquisition  system.  The  CCD  camera  was  an  ImagerPro  X  manufac¬ 
tured  by  LaVision  Co  (http://www.lavision.de).  It  has  a  maximum  resolution  of  1600  x  1200  pixels 
and  a  14-bit  pixel  dynamic  range.  The  maximum  frame  rate  at  full  resolution  is  15  Hz.  Image  data 
streams  into  a  dedicated  computer  system  through  a  CameraLink  interface,  and  is  directly  saved  to  a 
computer  hard  drive.  In  this  study,  a  frame  rate  of  10  Hz  was  selected,  along  with  an  image  resolu¬ 
tion  of  1200  X  1200  pixels.  The  Q-switched  laser  was  a  Solo  XT  PIV  made  by  New  Wave  Research 
(http://www.new-wave.com).  It  is  a  laser  diode-pumped  Nd:Yag  green  laser  (light  wavelength  =  532 
nm)  with  a  maximum  repetition  rate  of  15  Hz  (15  pairs  of  pulses  per  second).  The  duration  of  each 
laser  pulse  is  20  ns  and  the  light  energy  is  200  mJ  per  pulse.  A  programmable  timing  unit  (PTU)  was 
used  to  send  controlling  signals  to  trigger  and  synchronize  the  laser  and  camera.  Image  acquisition 
and  timing  controls  are  integrated  into  the  LaVision  software,  DaVis  7.2. 

Fluorescein  was  selected  as  a  surrogate  for  MC  because  fluorescein:  I)  exists  in  both  solid  and 
aqueous  phases  under  room  temperature,  and  2)  is  measurable  by  a  fluorometer  with  high  accuracy. 
Subsequently,  the  study  shows  that  the  molecular  diffusivity  (kinematic)  of  fluorescein  in  water  at 
20”  C  is  5.05x10  ''^  (m^/s)  is  on  the  same  order  of  magnitude  as  that  for  TNT  (lO  '**  m^/s),  both  being 
much  smaller  when  compared  with  the  hydrodynamic  diffusion  (~10'^  m^/s)  for  the  smallest  breach 
size  and  current  investigated  of  0.05  m/s,  which  governs  the  release  rate.  Therefore,  for  the  cases  we 
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studied,  the  molecular  diffusivities  of  MC  are  not  directly  relevant  to  the  study  results  but  are  still 
considered  appropriate  for  modeling  the  release  rates  of  MC. 

The  2-D  instantaneous  concentration  distribution  of  fluorescein  solution  was  measured  using 
planar  laser  induced  fluorescence  (LIF).  The  molecular  diffusivity  (kinematic)  of  fluorescein  in  water 
at  20°C  is  5.05X  10''®  (m^/s).  The  peak  absorption  of  fluorescein  occurs  at  494  nm  and  the  emission 
peak  at  521  nm.  A  DPSS  CW  laser  (www.laserglow.com)  with  wavelength  of  473  nm  and  300  mW 
of  light  power  was  applied  to  excite  the  dye  solution.  The  laser  beam  was  steered  by  an  optical 
scanner  (Cambridge  Technology  621  OH  galvanometer)  to  cut  across  the  dispersing  fluorescein  plume 
emitted  from  the  breach.  The  resulted  laser  “sheet”  was  perpendicular  to  the  direction  of  the  mean 
current  in  the  open  channel,  and  placed  at  a  location  downstream  of  the  breach  hole  (20-30  cm) 
where  the  plume  was  picked  up  by  turbulent  eddies,  grown  to  a  larger  size  and  sufficiently  diluted. 
This  was  considered  an  important  step  to  increase  the  accuracy  of  concentration  measurement,  as  the 
linear  relation  between  the  fluorescein  concentration  and  the  emitted  light  intensity  is  optimal  at 
lower  concentrations.  In  addition,  laser  light  attenuation  becomes  less  of  an  issue  at  lower 
concentrations.  A  mirror  was  placed  downstream  of  the  laser  sheet,  forming  a  45 -degree  angle  from 
the  mean  flow  direction.  A  CCD  (Point  Grey  Research,  FI)  IEEE  1394  camera  with  1032  x  776  pixel 
resolution  was  introduced  from  the  channel  sidewall  to  collect  images  of  the  fluorescent  light  emitted 
from  the  dispersing  dye.  An  USB-based  analog  output  module  (Data  Translation  DT9810)  was  used 
to  generate  signal  waveforms  to  synchronize  the  shutter  of  the  CCD  camera  and  the  scanning  angle 
of  the  galvanometer. 

3.3.3  Procedures 

The  first  set  of  experimental  studies  was  developed  to  examine  the  dominant  factors  that  determine 
the  mass  release  rate  from  a  breached  hole.  According  to  our  hypothetical  analysis,  the  three 
controlling  parameters  were  the  radius  of  the  semispherical  cavity  (R),  the  radius  of  the  crack  hole 
(b),  and  the  shear  velocity  above  the  breached  surface  (u*),  where  the  shear  velocity  over  a  smooth 
wall  has  a  well-known  relation  with  the  ambient  flow  speed  (U).  We  selected  three  cavity  radii,  R  = 
2.0,  3.5,  and  5.8  (cm)  and  five  hole  radii,  b  =  0.5,  2.5,  5.0,  10.0,  and  25.0  (mm),  respecttively. 
Repeatable  flow  conditions  were  established  with  three  nominal  free  stream  speed,  i.e.,  U=l,  15, 
and  3 1  cm/s,  and  the  corresponding  shear  velocities  were  about  {u*  =  3.4,  6.7,  and  1 1 .2  mm/s).  These 
were  referred  to  as  ’’Low  Current,”  “Median  Current,”  and  “High  Current”  conditions  hereafter.  As  a 
first-order  approximation,  we  hypothesized  that  the  time  scale  of  mass  depletion  was  related  to  the 
volume  of  the  cavity  (F=  2l3nl^),  the  area  of  the  hole  {A  =  nb^),  and  the  shear  velocity  (m*).  For 
dissolved  MC  (no  solid  MC  replenishment),  a  characteristic  time  to  depletion  can  be  defined  as 


V  _  2R^ 

AUf  3utb^ 


(14) 


In  this  study,  the  combination  of  the  selected  R  and  b  forms  a  wide  range  of  characteristic 
depletion  times  (from  8  to  77,657  seconds).  Table  1  shows  the  series  of  R  and  b  and  the 
corresponding  to  under  the  shear  velocity  u*  =  6.7  mm/s  (the  median  current  condition). 
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Table  1.  Characteristic  depletion  time  (seconds)  under  the  combination  of  cavity  radius  and  hole 
radius  under  the  median  current  condition. 


to 

b  =  0.5 
(mm) 

b  =  2.5 
(mm) 

b  =  5  (mm) 

b  =  10  (mm) 

b  =  25  (mm) 

R  =  2.0  cm 

3,184 

127 

32 

8 

N/A 

R  =  3.5  cm 

17,065 

683 

171 

43 

7 

R  =  5.8  cm 

77,657 

3,106 

777 

194 

31 

The  order  of  measurements  with  different  combination  of  b  and  R  were  randomized.  First,  the 
channel  flow  was  set  up  at  a  desired  current  speed  that  was  kept  unchanged.  PIV  was  then  applied  to 
characterize  both  the  near-breach  external  flow  and  the  induced  in-cavity  circulation  (Figure  12  and 
Figure  13).  A  pulsed  laser  sheet  was  then  delivered  from  the  top  of  the  channel.  The  light  sheet  was 
parallel  to  the  mean  flow  direction,  perpendicular  to  the  breached  surface,  and  cut  through  the 
semispherical  cavity  in  the  middle.  Since  the  physical  model  is  made  of  optically  clear  acrylic,  it 
allowed  excellent  optical  viewing  (for  both  the  laser  light  and  the  imaging  device)  to  the  inner  space 
of  the  cavity.  Profiles  of  the  external  mean  current  speed  U{z)  and  the  Reynolds  shear  stress  -m'w'(z) 
were  determined  from  ensemble  averaging  instantaneous  velocity  fields  measured  more  than  10 
minutes.  Thus,  the  shear  velocity  u*  could  be  estimated  by  either  fitting  the  profile  of  U  with  the  log- 
law  (hydraulically  smooth  boundary),  or  taking  the  square  root  of  the  maximum  value  of  the  -u'w' 
profile.  Measurements  on  internal  circulations  were  followed  to  determine  the  temporal  and  spatial 
variation  of  circulation  patterns  and  magnitudes. 
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Figure  12.  PIV  measurement  (with  laser  sheet). 
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Figure  13.  Measured  internal  flow  circulation  from  the  PIV  images. 


After  the  PIV  measurements,  the  cavity  model  was  taken  out  and  refilled  with  a  well-mixed 
fluorescein  solution.  Then  the  crack  hole  was  sealed  with  waterproof  tape  and  the  cavity  was 
mounted  back  on  to  the  platform.  LIT  image  sequence  acquisition  was  then  initiated  quickly  after  the 
seal  was  removed  and  a  fluorescein  plume  started  to  emit  from  the  breach  opening.  The  instantan¬ 
eous  concentration  distribution  on  a  2-D  cross-flow  downstream  section  was  measured  using  LIF  at  a 
rate  of  1  Hz  over  400  seconds  (400  images)  (Figure  14).  In  essence,  the  time  scale  of  mass  depletion 
is  tabulated  in  Table  1.  Depending  on  whether  the  depletion  time  was  much  smaller  than,  greater 
than,  or  the  same  order  as  the  duration  of  LIF  measurements,  the  mass  depletion  rates  are  estimated 
with  different  methods,  which  will  be  explained  in  Section  3.4.4. 

Calibration  measurements  were  also  performed  to  relate  the  measured  concentration  to  that  of  the 
source,  under  the  assumption  that  the  source  was  initially  well  mixed  before  it  was  injected  in  the 
modeled  cavity.  Using  a  pipette  (5  ml),  samples  of  dye  solution  were  withdrawn  from  the  cavity 
before  and  after  each  LIF  measurement.  The  samples  were  then  injected  into  a  small  aquarium  and 
diluted  1700  times  by  adding  clean  water  into  the  aquarium.  The  concentration  of  the  diluted  source 
was  then  measured  by  LIF  using  the  same  optical  configuration  as  in  flume  experiments.  This 
concentration,  modified  by  the  dilution  factor  (1700x)  is  the  concentration  of  the  source,  which  was 
used  to  normalize  the  plume  concentration  measurement  and  to  estimate  the  total  mass  depleted 
through  the  hole  during  the  flume  LIF  measurement. 
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Figure  14.  Measured  fluorescein  concentration  processes  from  the  LIF  images  on  a  cross-section 
downstream  across  the  plume. 

3.3.4  Estimation  of  the  Mass  Release  Rate,  F 

The  instantaneous  chemical  mass  release  rate  can  be  experimentally  determined  by  measuring  the 
streamwise  release  of  dye  concentration  at  a  downstream  cross-section.  In  theory,  the  mass  release 
rate  is 


^  =  J  J  ^  u{y,z)c{y,z)dydz ,  (15) 

where  A  is  the  downstream  cross  section,  y  and  z  are  the  spanwise  and  vertical  coordinates,  and  u 
and  c  are  fluctuating  streamwise  velocity  and  dye  concentrations.  According  to  Reynolds 
decomposition,  they  can  be  written  as 

u{y,z)  =  U(z)  +  u'{y,z) 

c(y,z)  =  C(y,z)  +  c'(y,z),  ^  ^ 

where  the  U  is  the  ensemble  mean  velocity,  which  is  assumed  to  be  homogenous  along  the 
spanwise  direction,  thus  only  a  function  of  z,  and  u  ’  is  the  turbulent  fluctuation.  Similarly,  C  is 
the  ensemble  mean  concentration  and  c  ’  is  the  fluctuation.  Thus,  it  can  be  shown  that  for  a 
statistically  steady  state,  the  temporally  averaged  total  flux  rate  is 


j  [^(^)(|  C(y,z)Jy)J/z+  Ij ^u'{y,zy{y,z)dydz, 


(17) 


where  the  two  integrals  on  the  right-hand  side  are  the  mean  current  flux  and  the  turbulent  flux, 
respectively.  Previous  studies  in  the  flume  laboratory  indicate  that  the  turbulent  flux  in  the 
streamwise  direction  is  generally  an  order  of  magnitude  smaller  than  the  mean  flux.  Therefore, 
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we  assume  the  turbulent  flux  in  our  measurement  to  be  »  mean  current  flux,  and  use  the  mean  flux 
to  estimate  the  mass  release  rate.  Note  that  it  is  technically  impractical  to  measure  the  streamwise 
turbulent  flux  over  a  cross-stream  section,  although  it  can  be  done  over  a  plane  that  is  parallel  to  the 
mean  flow. 

For  the  MC  release  rate  measurements  reported  here,  the  mass  release  rate  is  not  statistically 
steady,  especially  for  the  characteristically  small  depletion  time  tD.  The  temporal  variation  of 
releaserate  cannot  be  avoided  in  these  cases,  as  there  is  no  dissolution  interface  in  the  currently 
designed  system.  Based  on  the  argument  above,  the  instantaneous  release  rate  was  approximated  by 
ignoring  the  correlation  between  the  fluctuating  part  of  the  velocity  and  the  instantaneous 
concentration,  so  we  have 


^(0  =  \  '^^]u{z)c{y,z)dydz,  (18) 

where  H  is  the  height  of  the  image  above  the  flume  bed  and  W  is  the  width  of  the  image. 

The  mass  release  rate  was  determined  as  statistically  steady  over  the  period  of  image  acquisition 
(i.e.,  without  obvious  longterm  trend)  if  the  to  is  large  {e.g.,  >  400  s).  Otherwise,  the  observed  mass 
release  decreases  in  a  log-linear  fashion  with  time.  Figure  15  shows  a  time  series  of  F(t)  for  typical 
small,  medium,  and  large  values  of  fe.  Note;  for  small  and  medium  fc,  the  figure  is  presented  in  a 
semi-log  fashion.  Dashed  lines  are  fitted  to  exponential  functions  with  respect  to  time. 


u.  =  1 .2  (cm/s),  R  =  35  (mm),  b  =  25.0  (mm),  Iq  =  4  (s) 
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u,  =  1.1  (cm/s),  R  =  35  (mm),  b  =  2.5  (mm),t0  =  417  (s) 


Figure  15.  Time  series  of  measure  mass  release  rate  with  small,  median  and  large  depletion  time 
scale,  to- 
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For  scenarios  with  large  values  of  to,  like  the  third  scenario  in  Figure  15,  we  can  assume  the  mass 
release  rate  is  quasi-steady  state.  Thus,  the  ensemble  average  of  the  recorded  series  of  F{f)  is 
considered  as  the  MC  release  rate,  and  we  can  assume  a  model  of  the  mass  release  rate  as 

F  =  UfAC^,  (19) 

where  A  =  nb^,  the  hole  area,  Co  is  the  average  concentration  within  the  cavity,  and  w/is  the  “release 
velocity.”  The  release  velocity  has  previously  been  determined  empirically  as  0.023  C  according  to 
FLUENT  simulations  (with  Sc=  1),  and  0.025  t//5'c  according  to  lab  experiments  with  small  hole 
radius.  Our  new  data  show  that  m/ could  be  a  function  of  the  hole  radius  and  the  cavity  radius  as  well. 

For  smaller  values  of  to,  the  observed  log-linear  trend  is  characteristic  of  a  simple  first  order  mass 
depletion  model,  which  can  be  applied  to  describe  the  MC  release; 

r^  =  -F^-u,AC„  (20) 

where  Co  is  the  average  concentration  in  the  cavity,  as  it  decreases  with  time.  Eq.  (20)  indicates  that 
the  decreasing  release  rate  can  be  modeled  as 


F  =  Uj^ACq  exp 

V 


(21) 


where  Cq  is  the  initial  in-cavity  concentration,  which  can  be  determined  through  the  calibration 
experiment.  Linear  regression  can  be  applied  to  log  (F)  vs.  t  to  determine  the  release  speed  Uf. 

3.3.5  Shear  velocity  ih 

We  have  argued  that  mass  release  rate  from  a  small  breach  should  be  scaled  with  the  surrounding 
flow  structure  that  has  a  comparable  scale  as  the  breach  hole  itself  Therefore,  the  shear  velocity 
above  the  breached  surface  is  selected  as  one  of  scaling  parameters.  A  relationship  between  the 
ambient  current  speed  and  the  shear  velocity  can  be  established  based  on  the  roughness  and 
orientation  of  the  surface.  Therefore,  it  will  be  convenient  to  separate  the  Reynolds  number  of  the 
overall  flow  and  other  factors,  so  that  we  can  focus  on  the  fundamental  mechanisms  that  control  mass 
release  through  a  small  opening. 

In  our  previous  studies,  mass  release  rate  has  been  associated  with  the  depth-averaged  velocity  (or 
the  ambient  current  speed).  Here  we  present  the  relationship  between  shear  velocity  and  ambient 
speed  under  the  scenario  of  simple  uniform  flows  over  a  smooth  surface.  Most  of  our  experiments 
were  conducted  in  three  groups.  For  each  group,  three  nominal  flow  speeds  were  established:  low, 
median,  and  high.  Vertical  profiles  of  mean  velocity  and  Reynolds  shear  stresses  of  the  nine  flows 
are  shown  in  Figure  16. 

Shear  velocity  u*  was  estimated  as  the  square  root  of  the  peak  value  of  -u'w' .  The  ambient  current 
speed  U  was  estimated  as  the  depth-average  of  the  vertical  profile  of  the  mean  velocity.  These  values 
and  the  ratio  u*IU  are  listed  in  Table  2. 

Our  previous  study  has  indicated  that  for  general  cases,  the  mass  release  velocity  is  Uf=  0.025  U. 
With  the  new  scaling  parameter  u*,  we  should  expect  that  Uf=  0.50-0. 63m*,  since  u*IU  is  about  0.04 
to  -0.05  (Table  3). 
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3.3.6  Mass  Release  Velocity  Uf  vs.  Shear  Velocity  ih 

For  a  large  depletion  time,  measured  mass  release  rate  F  was  normalized  by  the  calibrated  source 
concentration  and  divided  by  the  crack  area  to  deduce  the  mass  release  velocity  Uf.  For  a  smaller 
depletion  time,  first-order  depletion  model,  log-linear  regression  analysis  based  on  Eq.  (20)  was 
directly  applied  to  deduce  the  release  speed  Uf.  For  the  median  depletion  rate,  both  methods  were 
applied  and  the  average  of  the  two  results  were  chosen  as  the  mass  release  velocity. 

First,  the  mass  release  velocity  was  plotted  against  the  shear  velocity  for  all  combinations  of  R  and 
b,  as  shown  in  Figure  17.  Generally,  Wf  increases  almost  linearly  with  u*  for  most  cases.  If  we  apply 
linear  regression  of  Uf~u*  for  each  pair  of  R  and  b,  and  force  the  intercept  to  be  zero,  the  scaling 
between  U{  and  u*  can  be  determined,  and  the  results  are  listed  in  Table  4. 


Figure  16.  Vertical  profiles  of  the  mean  velocity  and  the  Reynolds  shear  stress. 


Table  2.  Shear  velocities  and  ambient  current  speed  (cm/s). 


u* 

Low  Speed 

U 

u,/U 

u* 

Median  Speed 

U 

uJU 

u* 

High  Speed 

U 

u^/U 

Group  1 

0.328 

6.756 

0.049 

0.667 

15.692 

0.043 

1.152 

31.306 

0.037 

Group  2 

0.351 

7.113 

0.049 

0.663 

15.943 

0.042 

1.097 

27.087 

0.041 

Group  3 

0.330 

7.133 

0.046 

0.649 

15.360 

0.042 

1.111 

30.290 

omi 

Table  3.  Scaling  between  the  mass  release  speed  and  the  shear  velocity,  Uf /u*. 


Uf/U* 

b  =  0.5  (mm) 

b  =  2.5  (mm) 

b  =  5  (mm) 

h  =  10  (mm) 

b  =  25  (mm) 

R  =  2.0  cm 

0.44 

0.49 

0.51 

0.50 

N/A 

R  =  3.5  cm 

0.84 

0.90 

0.97 

1.25 

1.17 

R  =  5.8  cm 

1.68 

1.69 

1.62 

1.97 

1.69 

29 


Figure  17.  The  relation  between  mass  release  speed  Uf,  and  shear  velocity  u*. 

According  to  our  empirical  results,  Uf  =  0.44-1  .97m*,  a  much  wider  variability  than  our  preliminary 
simulation  results. 

3.3.7  Effect  of  Crack  Hole  Radius  b  on  the  Release  Speed  Uf 

According  to  the  values  listed  in  Tables  4  and  5,  the  scaling  factor  between  uf  and  u*  seems  to 
slightly  increase  with  the  crack  hole  radius  for  all  the  different  cavity  radii  R.  By  grouping  measured 
uf/u*  with  the  hole  radius  b,  a  single  factor  ANOVA  shows  that  the  dependency  of  uf/u*  on  the  hole 
radius  is  not  significant  (P -value  =  0.93)  (Table  4). 


Table  4.  ANOVA  results  for  Uf/u*  dependence  on  hole  radius. 


SUMMARY 


Groups 

Count 

Sum 

Average 

Variance 

B  =  0.5 

3 

2.96 

0.986666667 

0.400533333 

b  =  2.5 

3 

3.08 

1 .026666667 

0.372033333 

B  =  5 

3 

3.1 

1 .033333333 

0.311033333 

B  =  10 

3 

3.71 

1 .236666667 

0.540233333 

B  =  25 

2 

2.81 

1.405 

0.11045 

ANOVA 

Source  of 
Variation 

SS 

df 

MS 

F 

P-value 

Fcrit 

Between 

Groups 

0.305254762 

4 

0.07631369 

0.20452631 1 

0.929521104 

3.633088512 

Within  Groups 

3.358116667 

9 

0.373124074 

Total 

3.663371429 

13 
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3.3.8  Effect  of  Cavity  Radius  /?on  the  Release  Speed  Uf 

According  to  the  values  listed  in  Tables  4  and  5,  the  eavity  radius  has  a  strong  effect  on  the  sealing 
factor  between  Mf  and  u*.  By  grouping  measured  Uflu*  ratio  with  eavity  radius  R,  a  single  factor 
ANOVA  verifies  that  the  non-dimensional  u^u*  increases  with  R  significantly  (P-value  =  1.28  x  10'^) 
(Table  5). 

3.3.9  Observed  Universal  Function  for  Mass  Release  Speed 

If  we  exclude  the  hole  radius  b,  a  simple  dimensional  analysis  shows  that  the  normalized  mass 
release  speed  Uf/u*  should  be  a  function  of  the  cavity  radius  based  Reynolds  number: 


Re^ 


u»R 

? 

V 


(22) 


where  vis  the  kinematic  viscosity.  Figure  18  shows  the  scaling  factor  Uf/u*  versus  ReR  for  all  the 
measurements.  A  linear  regression  analysis  indicates  that 


-^  =  0.0042  Re^ -1-0.0492,  (23) 

with  the  correlation  coefficient  R^  =  0.7227. 

Table  5.  ANOVA  results  for  Uf  /  u*  dependence  on  cavity  radius. 


SUMMARY 


Groups 

Count 

Sum 

Average 

Variance 

R  =  10 

4 

1.94 

0.485 

0.000966667 

R  =  35 

5 

5.12 

1.024 

0.03003 

R  =  50 

5 

8.6 

1.72 

0.02035 

ANOVA 


Source  of 
Variation 

ss 

df 

MS 

F 

P-vaiue 

F  crit 

Between 

Groups 

3.458951429 

2 

1.729475714 

93.06444016 

1. 27801 E-07 

3.982297957 

Within  Groups 

0.20442 

11 

0.018583636 

Total 

3.663371429 

13 

3.4  EMPIRICAL  STUDY  -  BURIED  SHELL 
3.4.1  Pore-water  Flow  Experiment 

A  pore-water  tank  was  designed  and  fabricated  to  study  the  release  rate  of  tracer  dye  from  a  buried 
breached  dye  container  to  simulate  MC  from  buried  shells.  The  tank  includes  two  constant  depth 
water  reservoirs  and  a  connecting  porous  media  channel.  The  overall  configuration  and  dimension  of 
the  tank  is  shown  in  Figure  19. 

Hydraulie  pressure  drove  the  seepage  flow  through  the  porous  media  ehannel  and  was  aehieved  by 
maintaining  eonstant  water  levels  in  the  two  reservoirs.  Water  was  eontinuously  pumped  into  the 
high-head  reservoir  and  drained  through  one  of  the  three  holes  to  maintain  a  desired  water  level.  As 
water  seeped  through  the  porous  media,  it  upwelled  into  the  low-head  reservoir  and  eventually 
drained  through  the  hole  (see  Figure  19(b)).  The  system  was  designed  such  that  head  difference  on 
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the  two  sides  of  the  porous  media  channel  can  be  1,  2,  or  3  inches,  respectively.  Throughout  this 
study,  the  resultant  seepage  flow  conditions  are  called:  “low,”  “medium,”  and  “high”  flow 
conditions,  respectively. 

Figure  20  shows  the  operation  of  the  pore-water  tank.  The  head  difference  used  in  this  case  is 
3  inches.  Two  types  of  clean  sand  were  used  to  model  the  coastal  water  sediment:  fine  sand 
(mean  diameter  ds=  150  pm)  and  coarse  sand  (mean  diameter  Js  =  500  pm).  The  estimated 
parameters  for  the  two  types  of  sand  are  given  in  Table  6  (the  estimation  of  conductivity  will  be 
detailed  in  the  following  paragraphs). 


Figure  18.  Normalized  release  speed  vs.  cavity  radius  based  Reynolds  number. 


(a) 


Constant  water  head  kept  between  the  two  ends  to 
drive  the  pore  water  flow 


Figure  19.  (a)  Schematic  of  the  pore-water  tank;  (b)  Drain  holes  on  the  sidewalls  of  the  high-head 
and  low-head  reservoir. 
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Table  6.  Physical  properties  of  sand  in  the  dye  release  experiments. 


Parameter 

Fine  sand 

Coarse  sand 

(150  pm) 

(500  pm) 

Porosity 

31% 

28% 

Hydraulic  conductivity  (cm/s) 

0.011 

0.065 

Figure  20.  Dye  tracing  in  the  pore-water  flow  tank. 


Pore- water  seepage  velocity  was  measured  through  dye  visualization  (Figure  20),  and  was  found  to 
be  slightly  non-uniform  over  the  depth  of  the  channel.  Under  a  typical  medium  flow  condition,  the 
seepage  velocity  was  about  Us  =  0.20  mm/s  near  the  top  wall,  and  decreased  with  depth  until  1.2 
inches  below  the  top  wall.  It  then  became  nearly  uniform  (v  =  0.16  mm/s)  over  a  1.5-inch  thickness 
above  the  bottom  wall.  The  vertical  profile  of  the  seepage  velocity  was  similar  under  all  flow 
conditions.  During  the  dye  release  experiments  described  in  the  next  section,  the  modeled  shell  was 
buried  at  ~1  inch  above  the  bottom  wall.  The  local  seepage  velocity  was  then  estimated  using  the 
measured  total  discharge  and  the  estimated  vertical  velocity  profile, 

M,  =  0.88Q/(^f),  (24) 

where  the  Q  is  the  total  discharge,  fis  the  porosity,  A  the  area  of  the  channel  cross-section,  and  0.88 
was  the  mean  correction  coefficient  that  accounts  for  non-uniformity  of  all  the  flow  profiles. 

A  2-D  numerical  model  was  developed  to  simulate  the  flow  conditions  using  the  Matlab  Partial 
Differential  Equation  Toolbox  (see  Section  3.8  for  details).  The  model  solves  for  the  distribution  of 
piezometric  (hydraulic)  head  h  over  the  sandy  matrix.  The  Darcy  velocity  is  then  evaluated  as 
q  =  -KVh .  Figure  21  shows  the  simulated  result,  the  simulation  corresponds  closely  to  the 
experimental  observation  in  terms  of  the  visualized  flow  streamline.  It  also  shows  that  the  Darcy 
velocity  vectors  are  nearly  horizontal  in  most  regions  of  the  test  section. 
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Figure  22  shows  the  plot  of  the  averaged  head  gradient  function,  G  =  -W  f  {dhldx^iy  (where  if  = 

Jo 

3  in,  the  height  of  the  channel,  JV^  5  inches,  the  width  of  the  channel),  as  a  function  of  the  distance  x 
(from  0  to  14.5  inches).  It  reaches  a  nearly  constant  value  in  the  central  part  of  the  test  section, 
indicating  a  uniform  flow  condition.  Thus,  we  can  use  the  measured  discharge  Q  and  the  simulated 
result  G  (at  the  constant  value  of  the  horizontal  flow  test  section)  to  estimate  the  hydraulic 
conductivity,  AT  (  =  Q/G).  Using  this  approach  for  both  types  of  sediment,  we  found  =  0.01 1  (cm/s) 
for  the  fine  sand  and  K  =  0.065  (cm/s)  for  the  coarse  sand  (see  Table  6,  above). 


0  2  4  6  8  10  12  14 


X  (inch) 

Figure  21.  Numerical  model  simulation  result  of  the  flow  in  the  pore-water  tank. 


25 

20 

c<P  15 
E 

CD  10 

5 


X  (inch) 

Figure  22.  Averaged  head  gradient  function  plot. 

Table  7  lists  the  typical  discharge  values  for  porous  media  flow,  in  addition  to  the  estimated 
seepage  velocity  near  the  buried  shell  under  three  typical  flow  conditions  for  fine  and  coarse  sand 
structures. 
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Table  7.  Pore-water  flow  conditions. 


Flow 

conditions 

Head 

difference 

(cm) 

Discharge 

(ml/hour) 

Average 

Darcy 

velocity 

(v-mm/s) 

Average 
Seepage  Flow 
velocity 
(Usmm/s) 

Seepage  Flow 
velocity 
around  the 
shell 

(Us-mm/s) 

Fine 

Low 

2.54 

263 

0.008 

0.026 

0.023 

sand 

Medium 

5.08 

514 

0.015 

0.048 

0.042 

High 

7.62 

815 

0.023 

0.074 

0.065 

Coarse 

Low 

2.54 

1600 

0.046 

0.163 

0.143 

sand 

Medium 

5.08 

3250 

0.093 

0.332 

0.292 

High 

7.62 

4380 

0.125 

0.446 

0.393 

3.4.2  Dye  Release  Experiments 

Based  on  the  sandbox  porous  flow  (Figure  19  and  Figure  20),  the  dye  release  experiments  were 
conducted  as  visualized  in  Figure  23.  A  small  glass  bottle  (4-cm  high  and  2.5-cm  in  diameter,  inner 
volume  =11  ml)  with  a  small  circular  hole  drilled  in  the  cap  that  was  used  as  a  surrogate  for  a 
breached  shell  (Figure  24).  The  experimental  methodology  used  is  described  as  follows.  The  bottle, 
filled  with  Rhodamine  WT  fluorescent  dye  solution  (molecular  diffusion  coefficient  in  water  D  =  1.5 
X  10  ''  cm^/s)  of  known  concentration  (usually  Co  =  400  ppm),  was  placed  in  the  lower  portion  of  the 
porous  media  channel,  with  its  hole  at  5  cm  above  the  bottom  and  the  breached  surface  parallel  to  the 
mean  flow  direction.  As  the  dye  was  transported  with  the  pore-water  movement,  it  upwelled  into  the 
low-head  reservoir  where  it  was  mixed  by  air-induced  turbulence  at  the  water  surface  generated  from 
a  small  fan.  Finally,  well-mixed  dye  was  drained  from  the  low-head  reservoir  where  it  was  collected 
by  a  flask  and  remixed.  A  fluorometer  (Turner  Designs  Cyclops)  probe,  submerged  in  the  flask,  was 
used  to  measure  and  record  the  concentration  of  the  exit  dye  continuously  at  a  frequency  of  one 
sample  every  10  minutes.  Thus,  the  flux  rate  (F)  of  dye  (as  MC  surrogate)  can  be  estimated  as 


F  =  QC,.,,  (25) 

where  Q  is  the  discharge  of  the  groundwater  flow  and  Cout  is  the  exit  dye  concentration  in  the 
flask. 


Figure  23.  Dye  release  and  flux  experiments. 
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Figure  24.  The  container  with  a  breach  hole  (db  =  0.56  cm) 
used  as  an  experimental  breached  munition  surrogate. 


Three  breach  hole  diameters  were  used  in  these  experiments:  =  0.21,  0.4,  and  0.56  cm.  The  dye 
release  experiment  consisted  of  1 8  runs,  with  three  different  breach  diameters,  three  flow  conditions, 
and  two  types  of  sand.  Each  test  run  lasted  for  at  least  24  hours,  or  until  50%  of  the  total  dye  mass 
was  released  from  the  breach  hole. 

3.4.3  Mass  Depletion  Analyses 

A  dimensional  analysis  allows  us  to  postulate  the  following  MC  release  rate  model. 


F  =  CqAuj  , 


(26) 


where  Co  is  the  mass  concentration  near  the  breach,  A  is  the  area  of  the  breach,  and  w/is  a  mass 
release  velocity,  measured  experimentally.  A  similar  model  was  used  to  estimate  the  rate  of 
release  from  a  breached  hole  in  flowing  water. 

Once  the  dye  is  released  through  the  crack  of  the  container,  subsequent  transport  is  determined  by 
pore-water  advection  and  dispersion  in  the  sediment  channel.  It  is  also  assumed  that  the  fluorescent 
dye  in  the  container  is  well  mixed,  i.e.,  C  =  Co  everywhere  in  the  container.  This  assumption  is 
reasonable  as  the  mass  flux  is  expected  to  be  low  due  to  the  low  seepage  velocity,  thus  the  time  scale 
for  diffusion  within  the  cavity  is  short  compared  with  the  timescale  for  the  overall  mass  depletion, 
allowing  the  release  rate  to  be  expressed  as  the  following  simplified  equation: 


(27) 


where  V  is  the  volume  of  the  container.  Therefore,  the  concentration  decreases  exponentially  as 
follows: 


Co(0  =  C,.exp  -—Uft  , 

V  f  J 


(28) 


where  C  is  the  initial  dye  concentration.  Meanwhile,  the  concentration  of  dye  exiting  the 
reservoir  is  measured  in  addition  to  the  flux  (Eqs.  (25)  and  (26)),  resulting  in  =  Auj^Cq  .  If 

a  simple  first-order  release  model  is  correct,  we  would  expect  to  see  the  concentration  of  the  exit 
dye  also  decrease  exponentially. 


36 


Therefore,  we  have, 


=--V  +  l“ 


AUf  ^ 

^C, 


Q  ‘ 


(29) 


Eq.  (28)  represents  the  analytical  solution  of  the  dye  concentration  in  the  container,  which  dilutes 
exponentially  over  time,  whereas,  Eq.  (29)  represents  the  analytical  solution  for  the  exit  dye 
concentration  which  behaves  similarly  (exponentially  decaying)  to  the  dye  concentration  in  the 
container.  In  addition,  Eq.  (29)  can  also  be  expressed  explicitly  as  a  function  of  flow  parameters,  U{, 
Q,  breach  area  (A),  and  volume  (V)  of  the  container.  The  analytical  solution  of  the  exit  dye 
concentration  (Eq.  (29))  was  compared  with  that  measured  from  the  dye  release  experiments  in 
Section  4.4.2. 

3.5  MODEL  VALIDATION 
3.5.1  Simulation  Parameters 

The  algebraic  analytical  function  for  MC  release  was  validated  by  the  numerical  simulation  results 
using  the  FLUENT  model.  FLUENT  is  a  general  computational  fluid  dynamic  model  widely  used  in 
automobile  and  combustion  industries.  The  FLUENT  model  includes  a  mesh  generating  software, 
GAMIT,  which  is  flexible  in  generating  model  configuration  and  meshes  for  scenarios  designed  by 
the  users.  For  this  numerical  modeling  study,  a  numerical  flume  channel  was  used  that  was  1  m  in 
length,  0.4  m  in  width,  and  0.2  m  in  depth.  A  circular  hole  (radius  =  b)  was  “drilled”  into  the  bottom 
of  the  channel.  The  hole  was  a  representative  of  the  breach  hole  in  the  shell,  and  it  connected  to  a 
semi-sphere  space  (radius  =  R)  beneath  the  flow  channel  as  seen  in  Figure  25.  The  flume  channel  and 
the  semi-sphere  were  filled  with  water.  The  inner  surface  of  the  semi-sphere  was  defined  as  the  MC 
solid-aqueous  interface,  through  which  the  MC  dissolution  occurs.  Solution  inside  the  semi-sphere 
cavity  was  released  through  the  hole  due  to  the  current  in  the  channel.  The  coordinate  system  was  set 
such  that  the  center  of  the  hole  was  the  origin,  positive  x  was  in  the  mean  flow  direction,  y  was  in  the 
cross-channel  direction,  and  z  was  positive  upward  (normal  to  flow). 


Figure  25.  FLUENT  model  geomefry  and  meshes  generated  by  GAMIT. 
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Model  meshes  with  variable  sizes  are  generated,  with  an  average  grid  size  of  10  mm.  Meshes  are 
refined  near  the  breach  hole  and  the  boundaries  of  the  channel  bottom.  Meshes  are  then  imported  into 
FLUENT  for  numerical  simulations. 

A  uniform  in- flow  profile  with  a  velocity  U  is  applied  to  the  upstream  cross-section  of  the  channel. 
Flow  is  simulated  as  a  laminar  type  with  non-slip  wall  boundary  condition.  At  the  solid-liquid  inter¬ 
face  inside  the  semi-sphere  (r  =  R,  where  r  is  the  radial  coordinate  originating  from  the  center  of  the 
hole  pointing  radially),  the  MC  concentration  is  assumed  to  reach  saturation  (C  =  Cs).  The  molecular 
diffusivity  of  MC  solution  is  assumed  to  be  Du  =  1 .0  x.  1  O'  m^/s,  /.  e. ,  the  Schmidt  number  of  the  MC 
solution.  Sc  =  vIDm  =  1,  where  v  is  the  viscosity  of  water. 

Steady-state  solutions  for  both  velocity  field  and  concentration  field  are  obtained  by  solving  the 
segregated  Navier-Stokes  equation  and  scalar  transport  equations  in  FLUENT.  Concentration  flux 
induced  by  channel  flow  is  evaluated  by  integrating  the  transport  term  L/C  over  a  selected  down¬ 
stream  cross-section.  Under  steady-state  conditions,  the  flux  should  be  independent  of  location 
downstream  of  the  breach  hole.  This  is  verified  by  calculating  fluxes  at  four  different  cross-sections 
in  FLUENT.  Figure  26  shows  a  typical  simulated  concentration  field  displayed  on  the  central  vertical 
plane  of  the  channel,  a  cross-section  25  cm  downstream  of  the  breach  hole,  the  plane  near  the  bottom 
of  the  channel,  and  a  central  vertical  plane  inside  the  semi-sphere.  A  saturation  constant  of  100  mg/L 
is  assumed  in  this  numerical  model  run.  As  is  shown  in  Figure  26,  solutions  inside  the  semi-spherical 
source  container  retain  concentrations  ~13  mg/L.  Concentrations  outside  the  hole  are  0.2  to  0.5  mg/L 
which  is  rapidly  reduced  downstream  (-0.002  mg/L,  at  a  location  25  cm  downstream). 


Figure  26.  Contours  of  concentration  on  the  central  vertical  plane  and  a  cross-section.  In  this  case, 
the  breach  hole  is  b  =  5  mm  in  radius  and  the  current  speed  is  U  =  6  cm/s. 
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3.6  MODELING  FATE  AND  TRANSPORT  OF  TNT,  RDX,  AND  HMX 

The  release  rate  function,  developed,  and  validated  by  FLUENT  previously,  provides  release  rates 
for  a  single  shell  under  various  hydrodynamic  and  shell  integrity  conditions.  Since  we  do  not  have 
adequate  information  about  the  quantity  and  magnitude  of  breached  shells  in  DoD  coastal  and 
estuarine  waters,  the  following  assumptions  are  made: 

1 .  The  release  rate  function,  developed  for  this  study,  can  be  used  to  predict  release  rate  for  any 
single  breached  shell. 

2.  For  multiple  breached  shells,  the  total  release  rates  can  be  obtained  by  summing  the  release 
rates  of  each  individual  beached  shell.  This  requires  an  estimate  or  assessment  of  the  number 
and  type  of  shells  at  a  given  site. 

In  the  following  example,  we  used  two  models  to  simulate  the  fate  and  transport  of  MCs  released 
from  breached  shells,  the  hydrodynamic  and  transport  model,  TRIM2D  (Cheng,  Casulli,  and  Gartner, 
1993;  Wang  et  al,  1998a  and  1998b),  and  the  USEPA’s  standard  water  quality  model,  WASP? 
(Wool  et  al,  2001;  Ambrose  et  al,  1992).  We  used  San  Diego  Bay  as  the  test  site,  since  the 
hydrodynamics  are  adequately  known  and  relevant  data  are  available  for  modeling  purposes.  For 
both  models,  we  choose  a  theoretical  munition  shell  with  the  following  dimensions  shown  in  Table  9. 


Table  8.  Sizes  and  dimensions  of  the  test  breached  shell. 


Length  (inches) 

Diameter  (inches) 

Hole  radius  (mm) 

Average  current  (cm/s) 

MC  Mass  (kg) 

24 

8 

20 

8 

8 

For  TRIM2D,  which  is  a  2-D  hydrodynamic  and  transport  model,  time-varying  current  speeds, 
predicted  by  the  model,  are  input  into  the  release  rate  function  for  the  breached  shell.  The  resulting 
release  rates  are  time  varying,  fluctuating  with  ambient  current  speed.  For  WASP,  a  fixed  averaged 
current  of  8  cm/s  is  chosen,  which  is  close  to  the  tidally  averaged  current  speed  in  San  Diego  Bay. 
Release  rate  for  WASP?  is  assumed  to  be  a  constant. 

3.6.1  TRIM2D  Model 

TRIM2D  is  a  depth-averaged  tidal  and  residual  circulation  model  with  a  finite-difference 
numerical  grid  and  scheme.  TRIM2D  has  been  previously  applied  to  simulate  time-varying  water 
surface  elevation,  averaged  water-column  currents,  and  associated  transport  of  contaminants  for 
several  estuarine  systems,  including  San  Francisco  Bay,  California;  Boston  Harbor,  Massachusetts; 
Charleston  Harbor,  South  Carolina;  and  Venice  Lagoon,  Italy  (Cheng,  Causulli,  and  Gartner,  1993; 
personal  communication).  Over  the  past  two  years,  SSC  Pacific  has  applied  TRIM  to  several  water- 
resource  and  water  quality  studies  for  San  Diego  (Chadwick  et  al,  1999;  Wang,  Sutton,  Richter,  and 
Chadwick,  2000;  Wang,  Chadwick,  Johnson,  and  Grovhoug,  2006). 

The  use  of  a  2-D  depth-averaged  model  such  as  TRIM  is  justified  because  field  data  show  that  the 
flow  in  San  Diego  Bay  exhibits  strong  uniformity  in  the  water  column  (Wang  et  al,  1998a).  Such 
uniformity  of  flow  in  the  water  column  results  from  the  fact  that  San  Diego  Bay  is  shallow  and  flow 
in  the  bay  is  primarily  driven  by  tides  from  the  Pacific  Ocean.  Other  assumptions  in  model 
formulation  included  incompressibility  of  water  and  the  Boussinesq  approximation,  which  simplifies 
effects  of  density  difference  by  only  considering  buoyancy  (associated  with  gravity)  and  neglecting 
other  high-order  terms  of  density  differences.  As  with  any  depth-averaged  model,  velocity  and 
density  are  implicitly  assumed  as  nearly  constant  over  the  water  column.  However,  horizontal  density 
gradients  are  treated  explicitly  in  the  momentum  equations. 
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We  approximated  bottom  shear  stress  using  a  Manning-Chezy  formulation  with  Manning's  n 
coefficient  assigned  as  a  function  of  water  depth: 

Tb  =  pul ,  (30) 

where  Xb  is  the  shear  stress  at  the  bottom,  p  and  u*  are  water  density  and  shear  velocity, 
respectively. 

The  shear  velocity  is  correlated  with  the  Manning-Chevy  relationship  as  follows: 


g\v\v 

c 

(31a) 

^1/3 

(31b) 

n 


where  g  is  gravity,  V  is  the  current  velocity  (simulated),  C  is  the  Chevy  coefficient,  and  n  is 
Manning’s  n-value. 

The  empirically  based  Manning-Chevy  formulation  provides  a  closure  to  the  bottom  boundary. 

The  only  empirical  parameter  needed  for  TRIM2D  is  Manning’s  n-value,  which  is  assumed  to  be  a 
function  of  water  depth. 

An  additional  transport  equation  is  used  for  each  chemical  species,  to  simulate  solute  transport. 
Solutes  are  assumed  to  be  dilute  (no  more  than  1000  parts  per  million),  not  affecting  water  density, 
and  thus  the  solute  transport  equations  are  uncoupled  from  hydrodynamics.  Furthermore,  the  trans¬ 
port  equation  is  solved  one  time-step  behind  the  continuity  and  momentum  equations,  effectively 
uncoupling  the  transport  equation  (Cheng,  Casulli,  and  Gartner,  1993).  This  approach  is  valid 
because  baroclinic  forcing  (i.e.,  density  flow)  changes  less  rapidly  than  barotropic  forcing  (gravity 
flow).  The  advective  term  in  the  transport  equation  is  processed  using  an  upstream  differencing 
scheme  that  ensures  conservation  of  mass  as  solutes  are  transported  by  advection.  Similar  numerical 
schemes  were  used  in  several  previous  studies,  in  which  details  about  the  formulation  and  previous 
applications  of  TRIM2D  can  be  found  (Wang,  Mill,  Martin,  and  Wool,  1997). 

TRIM2D  calculations  reported  here  were  compared  with  the  1983  tide  data  collected  by  the 
National  Oceanic  and  Atmospheric  Association  (NOAA)  at  the  Ballast  Point,  Downtown  San  Diego, 
and  South  Bay  stations  for  San  Diego  Bay.  We  used  tide  data  collected  at  the  Scripps  Piers  as 
boundary  conditions  for  the  model.  These  data,  along  with  the  Manning’s  n-coefficients,  were  fine- 
tuned  to  fit  the  model-simulated  tides  against  measured  values.  Once  calibrated,  we  validated  the 
model  by  independently  simulating  tides  and  currents,  and  then  compared  simulated  results  with 
measured  tides  and  currents.  Overall,  the  observed  differences  between  model  and  measurements 
were  less  than  3%  for  tidal  heights  and  less  than  15%  for  tidal  currents.  The  differences  in  tidal 
currents  between  the  model  and  measurements  may  be  attributed  to  multiple  sources,  including  errors 
associated  with  tidal  measurement  and  the  approximation  and  simulation  of  3-D  flows  by  the  2-D, 
vertically  averaged  model  (TRIM2D). 

From  a  rest  condition,  we  ran  TRIM  for  2  days  (48  hours)  to  ensure  that  a  steady-state  flow 
condition  was  reached  before  the  constant-load  effluent  entered  into  the  bay.  This  initial  2-day  run 
was  necessary  to  eliminate  transient  flows  that  should  not  exist  under  steady-state,  quasi-repetitive 
tidal  conditions.  We  used  a  6-minute  time-step  in  the  model  and  simulated  tidal  height  and  tidal 
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currents.  Simulation  for  1  year  (365  days)  was  executed  and  results  at  the  end  of  365  days  were 
analyzed.  The  Modeling  Results  sections  of  this  report  discuss  the  findings  from  this  modeling  study. 

3.6.2  TRIM2D  Model  Setup 

A  2-D  rectangular  numerical  grid  system  (of  dimensions  100  m  x  100  m)  was  used  for  TRIM. 
These  numerical  grids  cover  the  entire  San  Diego  Bay  and  portions  of  open  ocean  outside  the  mouth 
of  the  bay,  consisting  of  30,845  grid  cells  with  21,563  water  cells.  Hydrodynamics  of  the  bay  is 
driven  by  ocean  tides  only  and  measured  tides  are  prescribed  along  the  west,  south,  and  north 
boundaries  of  the  model  domain,  which  cover  the  entire  ocean  boundary. 

Two  simulation  scenarios  for  TRIM  were  evaluated.  The  first  scenario  explores  the  release  of  MCs 
from  a  shell  located  at  the  south  end  of  the  inner  bay.  Figure  27,  which  resembles  the  loading 
scenario  for  the  WASP7  model  (Box  #1).  The  second  simulation  scenario  assumes  that  MC  release 
occurs  at  the  San  Diego  Naval  Station  near  the  middle  of  the  bay.  For  both  simulation  scenarios, 
current  speeds  at  both  locations  simulated  by  TRIM2D  are  used  for  estimating  release  rate. 

Therefore,  the  release  rates  are  temporally  dependent  (time-varying). 
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Figure  27.  Location  (X)  of  MC  release  from  a  single  breached  shell,  inner  San  Diego  Bay  and  San 
Diego  Naval  Station,  respectively  (no  loading  scenario). 

3.6.3  WASP7  Model 

The  second  model  applied  for  this  analysis  is  the  toxic  module  (TOXI4)  of  the  U.S.  EPA  Water 
Analysis  Simulation  Program  Version  7  (WASP7)  (Wool  et  al,  2001;  Ambrose  et  al,  1993). 

WASP7  is  a  dynamic  compartment  (“box”)  modeling  program  for  aquatic  systems.  The  time-varying 
processes  of  advection,  dispersion,  point  and  nonpoint  mass  loading,  and  boundary  exchange  are 
represented  in  the  basic  program.  For  each  box,  WASP7  calculates  the  mass  of  contaminants  exiting 
the  box  by  physical  transport  and  chemical  or  biological  transformation,  or  accumulation  therein. 

The  WASP  models  have  been  applied  previously  to  address  fate  and  transport  in  estuaries, 
including  the  eutrophication  in  Tampa  Bay  (Wang,  Martin,  Wool,  and  Morrison,  1999)  and  a 
pesticide  spill  in  the  Sacramento  River  (Wang  et  al,  1997).  Model  accuracy  and  stability  were 
demonstrated  in  these  studies  and  we  have  selected  this  model  for  simulating  the  fate  and  transport  of 
MCs  leached  from  a  breached  shell  in  San  Diego  Bay. 
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The  WASP7  model  domain  includes  both  the  water  column  and  sediment  bed.  It  is  assumed  that 
the  water  column  (vertically  from  water  surface  to  the  sediment-water  interface)  is  fully  mixed. 
Below  the  water  column  is  a  vertically  stratified  sediment  column.  The  sediment  is  well-mixed 
within  each  layer,  but  segmented  vertically  into  a  well-mixed  active  surface  layer  and  a  deep 
sediment  layer.  The  deep  sediment  is  further  segmented  into  contaminated  and  clean  sediment 
regions.  Transport  pathways  incorporated  into  the  WASP7  model  include  volatilization,  burial, 
resuspension,  settling,  advection,  in  addition  to  pore-water  diffusion,  sorption,  and  decay.  Spatial 
variations  are  described  and  evaluated  by  the  use  of  multiple  boxes  in  the  horizontal  plane. 

Pollutants  in  the  water  column  can  be  transferred  to  the  atmosphere  by  volatilization  and  to  the 
sediment  layer  by  settling.  Concurrently,  pollutants  in  bottom  sediments  can  be  released  into  the 
water  column  by  resuspension  of  particles,  mixing  by  benthic  organisms,  and  diffusion  from  the 
sediment  pore  water. 

The  parameters  described  above  for  use  with  WASP7  are  listed  in  Table  9. 


Table  9.  Model  parameters  for  WASP7  (San  Diego  Bay  and  Elizabeth  River). 


Water  Column 

Values 

San  Diego  Bay 

Elizabeth  River 

Porosity 

0.8 

0.8 

Particle  specific  gravity 

2.65 

2.65 

Suspended  Solids 

Concentration 

2  mg/L 

15  mg/L 

Weight  Fraction  Carbon  in 

Solids 

0.02 

0.02 

Settling  velocity 

0.5  m/hr 

4.2  m/hr 

Resuspension  velocity 

9.4E-05  m/yr 

0.25  m/yr 

Molecular  diffusivity 

1.0.E-09  m^/s 

1 .0.E-09  m^/s 

Henry’s  Constant 

0.00013700  atm-m®/gmole 

0.00013700  atm-m^/gmole 

Molecular  weight  (TNT) 

227  g/mole 

227  g/mole 

These  parameters  are  generated  based  on  the  best  available  knowledge  that  we  have  for  San  Diego 
Bay  and  the  Elizabeth  River.  With  these  parameters,  the  fate  and  transport  of  MCs  from  a  single 
breached  shell  can  be  simulated.  Steady-state  solutions  were  obtained  from  the  WASP  model  runs. 

3.6.4  WASP  Model  Setup 

Two  model  scenarios  for  two  Navy  harbors  were  set  up  for  this  analysis.  The  first  model  scenario 
was  for  San  Diego  Bay,  California,  and  the  second  model  scenario  for  Elizabeth  River,  Virginia.  San 
Diego  Bay  is  an  estuary  in  which  the  flows  are  governed  by  tidal  forcing  from  the  Pacific  Ocean.  The 
bay  receives  very  little  freshwater  inflow  throughout  the  year  except  during  the  storm  seasons  of  the 
wet  years.  Transport  and  flushing  of  the  bay  are  primarily  driven  by  tides.  In  comparison,  the 
Elizabeth  River  watershed  is  a  riverine  sub-estuary  of  the  Chesapeake  Bay  estuarine  systems.  The 
watershed  is  a  complex  system  of  branches  and  small  tributaries  that  also  contains  a  deep-water 
harbor  used  by  the  U.S.  Navy  and  many  maritime  industries.  The  watershed  is  approximately  600 
km^,  with  ambient  water  and  sediment  quality  influences  from  the  highly  industrialized  commercial 
and  military  activities  and  densely  populated  urban  areas  of  the  metropolitan  Hampton  Roads 
Region.  The  river  flow  is  affected  by  the  tidal  motion  from  the  river  mouth  and  the  riverine  flows 
entering  from  upstream  of  five  major  watersheds.  As  a  result,  transport  and  flushing  of  the  river 
behave  differently  than  those  for  San  Diego  Bay. 
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San  Diego  Bay,  California 

An  eight-box  model  grid  for  San  Diego  Bay  was  generated,  as  shown  Figure  28.  Information  on 
the  grid  size  is  ineluded  in  Table  10  and  Table  11.  Flows  and  properties  in  eaeh  of  the  boxes  are 
assumed  to  be  homogeneous  (fully  mixed).  Transport  in  the  water  column  includes  net  advective 
flow  from  the  inner  bay  (Box  #1)  toward  the  mouth  of  the  bay  (Box  #8),  and  hydrodynamic 
dispersion  occurs  between  boxes.  Total  suspended  solids  (TSS)  concentration  in  the  water  column  is 
obtained  based  on  the  field  data.  The  settling  velocity  of  TSS  is  used  to  describe  sedimentation 
processes  (TSS  and  particulate  MCs).  The  resuspension  velocity  is  calculated  based  on  the  assigned 
settling  velocity  so  that  the  sediment  bed  remains  at  steady-state  condition,  i.e.,  there  is  no  gain  or 
loss  of  sediment  mass  in  the  bed.  The  hydrodynamic  dispersion  coefficients  between  adjacent  boxes, 
obtained  by  calibration  in  a  previous  study  (Wang,  Chadwick,  Johnson,  and  Grovhoug,  2006),  are 
assigned  to  the  two  model  grids. 


12  Kilometers 


Figure  28.  Model  grid  with  eight  boxes  for  WASP7  (box  numbering  starts  from  the  inner  bay  and 
ends  at  the  bay  mouth). 
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Table  10.  Grid  sizes  for  the  eight-box  model  grid  for  San  Diego  Bay. 


No. 

Surface  Area 
(m^) 

Perimeter  (m) 

Characteristic 
Length  (km) 

Water  Depth 
(m) 

Water  Volume 
(m^) 

1 

6,248.273.28 

12,733.26 

2.340 

2.0 

12,496,546.56 

2 

6,708,257.83 

15,821.94 

2.716 

3.3 

22,137,250.84 

3 

5,611,943.57 

10,562.35 

2.154 

4.8 

26,937,329.14 

4 

7,438,031.79 

13,831.00 

2.213 

6.0 

44,628,190.74 

5 

4,298,172.22 

11,254.83 

2.200 

6.6 

28,367,936.65 

6 

3,119,004.66 

12,988.63 

2.873 

9.2 

28,694,842.87 

7 

5,658,497.40 

18,893.49 

3.067 

9.8 

55,453,274.52 

8 

5,778,423.78 

18,386.98 

4.408 

10.2 

58,939,922.56 

Table  1 1 .  Simulation  matrix  for  San  Diego  Bay  using  WASP7. 


8-box 

San  Diego  Bay 

Hydrodynamic 

Dispersion 

Settiing/Resuspension 

Veiocities 

Mineralization 

Case  la 

Yes 

No 

No 

Case  1  b 

Yes 

No 

Yes 

Case  2a 

Yes 

Yes 

No 

Case  2b 

Yes 

Yes 

Yes 

To  understand  crucial  processes  governing  the  fate  and  transport  of  TNT  in  San  Diego  Bay,  we 
have  selected  the  following  four  cases,  each  with  different  combinations  of  model  parameters; 

Case  la  simulates  the  transport  of  TNT  in  San  Diego  Bay  driven  only  by  the  hydrodynamic 
dispersion.  This  scenario  is  similar  to  the  scenario  simulated  by  TRIM2D  with  no  mineralization. 

The  only  difference  is  that  San  Diego  Bay  is  represented  by  eight  boxes  in  WASP7  and  the  TRIM2D 
model  uses  much  finer  grid  resolutions.  With  the  different  resolutions,  WASP7  model  results  for 
Case  la  are  compared  with  those  from  TRIM2D;  both  can  be  considered  as  the  conservative 
estimates  of  TNT  in  water  column  for  this  analysis. 

Case  lb  is  the  same  as  Case  la,  but  with  mineralization  (same  as  the  rate  used  in  TRIM2D). 
Results  are  also  compared  with  the  results  from  the  corresponding  TRIM2D  run. 

Case  2a  simulates  both  transport  in  water  column  and  exchange  of  particulate  TNT  across  the 
water/sediment  interface  via  settling  and  resuspension.  Settling  velocity  of  0.5  m/hour  is  assumed, 
which,  according  to  the  Stoke’s  Law,  corresponds  to  particles  with  size  ~20  pm.  TSS  concentration 
of  2  mg/L  is  assumed  in  the  model  (derived  from  San  Diego  Bay  field  data).  Under  the  assumption 
that  sediment  mass  is  balanced  (no  net  sediment  erosion  or  accumulation),  a  resuspension  velocity 
necessary  to  satisfy  this  condition  is  calculated  and  assigned  to  the  model. 

Case  2b  is  the  same  as  Case  2a,  but  with  mineralization. 

For  Cases  2a,  2b  (and  later  in  the  4a  and  4b  cases  for  the  Elizabeth  River),  a  partitioning 
coefficient  between  the  particulate  and  aqueous  phases  of  TNT  is  needed.  A  literature  search  reveals 
little  or  no  direct  measurements  of  partitioning  coefficient  of  TNT  in  marine  water/sediment. 
Although  in  theory,  partitioning  coefficient  of  TNT  can  be  calculated  from  Kqw  and  TOC,  there  are 
studies  in  measuring  partitioning  coefficient  values  for  marine  water.  During  measurements  of  TNT 
partitioning  coefficients  in  laboratory  settings  by  Pennington  (personal  communication),  a  steady- 
state  distribution  between  sediment  and  overlying  water  could  never  be  reached.  “TNT  and  its 
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degradation  products  continued  to  decline  until  the  end  of  the  test.  This  can  be  attributed  to 
continuing  removal  of  the  compound  from  the  water  by  degradation,  transformation,  or  covalent 
bonding  to  sediment  components  such  as  organic  carbon”  (Brannon  et  al,  2005).  Mostly  likely,  these 
processes,  observed  by  Pennington  as  mentioned  above,  are  not  reversible,  contrary  to  partitioning 
that  is  used  in  the  model  as  a  reversible  process. 

In  other  efforts,  measured  partitioning  coefficients  of  TNT  in  soils  with  various  contents  of 
moisture  and  total  organic  carbon  were  reported  (Brannon  and  Pennington,  2002).  Gerald  and  Dortch 
(personal  communication)  summarized  and  examined  the  data  reported  by  Brannon  and  Pennington 
(2002)  and  confirmed  that  the  partitioning  coefficient  for  TNT,  RDX,  and  HMX  are  dependent  on  a 
number  of  soil  and  chemistry  parameters,  including  percentage  of  clay,  cation  exchange  capacity 
(CEC),  and  TOC.  Ranges  of  values  and  mean  of  the  partitioning  coefficients  associated  with  these 
parameters  are  shown  in  Table  12.  Measured  mean  TNT  partitioning  coefficients  (for  soil)  were  low, 
in  the  range  of -0.88  to  5.5  L/Kg.  Since  these  are  the  only  reported  values  that  are  conformal  to 
partitioning,  we  used  these  values,  0.9  and  5.0  L/Kg,  in  our  modeling  study,  despite  the  fact  that  they 
are  for  soil. 


Table  12.  Ranges  and  mean  values  of  soil  partitioning  coefficients  for  TNT  and  RDX  (Brannon  and 
Pennington,  2002). 


% 

Clay 

CEC, 

milliequivalents 

{g^>) 

TOC 

(%) 

Range  of 
partitioning 
coefficient  values 
L/Kg 

Mean  Partitioning 
Coefficient  Vaiues 

Number  of 
Observations 

TNT 

0-20 

0-10 

0-1 

1 .04-3.64 

0.88 

11 

20- 

50 

11-30 

1-3 

2.3-6.16 

3.39 

14 

>50 

>  30 

>  3 

2.23-11 

5.54 

12 

RDX 

0-20 

0-10 

0-1 

0.07-1.57 

0.53 

11 

20- 

50 

11-30 

1-3 

0.06-1.65 

0.85 

13 

>50 

>  30 

>  3 

0.31-8.4 

2.31 

15 

Elizabeth  River,  Virginia 

A  box  model  for  the  Elizabeth  River  was  generated,  as  shown  in  Figure  29.  The  box  model  is 
composed  of  26  boxes  for  the  water  column  and  the  sediment  layer,  boxes  1-13  and  boxes  14-26, 
respectively.  Box  #1  is  located  at  the  upstream  boundary  and  Box  #13  at  the  downstream  boundary 
for  the  Elizabeth  River.  Dimensions  of  the  model  boxes  are  shown  in  Table  13.  Flows  and  properties 
in  each  of  the  boxes  are  assumed  to  be  homogeneous  (fully  mixed).  Transport  in  the  water  column 
includes  net  advective  flow  from  the  upstream  segment  (Box  #1)  toward  the  mouth  of  the  river 
(Box  #13),  and  hydrodynamic  dispersion  between  boxes  (Figure  30).  Total  suspended  solids  (TSS) 
concentration  in  the  water  column  is  obtained  based  on  field  data  (Wang  et  al,  2007).  Settling  of 
particulate  MCs  are  described  by  the  settling  velocity  of  TSS.  The  resuspension  velocity  is  calculated 
based  on  the  assigned  settling  velocity  so  that  the  sediment  bed  deposition  rate  of  1  cm/year  is 
achieved,  which  is  obtained  from  sediment  coring  data  (Wang  et  al,  2007).  The  hydrodynamic 
dispersion  coefficients  between  adjacent  boxes,  obtained  by  calibration  in  a  previous  study  (Wang 
et  al.,  2006)  are  assigned  to  the  two  model  grids. 
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The  salinity  gradient  results  from  the  balance  of  advection  and  dispersion  in  the  river.  The 
freshwater  riverine  inflows  from  the  watersheds  are  obtained  from  measured  data  as  shown  in 
Table  14  and  Figure  31.  To  estimate  the  dispersion  coefficients,  we  used  the  finite-difference 
model  of  salinity  by  Thomann  and  Mueller  (1987)  to  estimate  the  dispersion  coefficients  for  the 
Elizabeth  River  (Figure  32). 


Box  Model  for  Elizabeth  River 


Segment  13 
Segment  12 
Segment  1 1 
Segment  10 
Segment  09 
Segment  08 
Segment  07 
Segment  06 
Segment  05 
Segment  04 
Segment  03 
Segment  02 
Segment  01 
Coastlines 


s 


Figure  29.  Box  model  for  Elizabeth  River. 
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Figure  30.  Conceptual  box  model  for  fate  and  transport  of  TNT  in  Elizabeth  River. 


Table  13.  Grid  sizes  for  the  13-box  model  grid  for  Elizabeth  River. 


Segment  No. 

Surface 

Area 

(m^) 

Depth 

(m) 

Volume 

(m^) 

Interaction  Area 
(m^) 

1 

261,214 

6.7 

1,750,410 

2,020 

2 

352,011 

5.0 

1,759,852 

2,768 

3 

359,596 

5.8 

2,081,851 

2,733 

4 

297,811 

7.3 

2,162,238 

2,392 

5 

366,058 

8.7 

3,198,533 

3,359 

6 

393,763 

6.8 

2,686,806 

4,465 

7 

444,106 

10.5 

4,658,982 

3,364 

8 

300,423 

8.0 

2,416,018 

3,801 

9 

248,628 

10.0 

2,497,280 

3,201 

10 

309,713 

7.5 

2,321,518 

2,898 

11 

529,690 

5.6 

2,992,211 

3,892 

12 

554,277 

6.8 

3,759,952 

3,497 

13 

654,112 

6.4 

4,188,029 

5,011 
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Table  14.  Discharge  data  for  Elizabeth  River. 


Watersheds 

Discharge 

(cms) 

Lafayette  River 

0.53 

Western  Branch 

0.94 

Eastern  Branch 

0.88 

Deep  Creek 

2.60 

Southern  Branch 

1.66 

Southern  Branch 
Avg.  Q=1.66cms 


Figure  31.  Riverine  freshwater  inflows  from  five  major  river  branches  of  the  Elizabeth  River 
watershed. 
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Figure  32.  Dispersion  coefficients  for  Elizabeth  River. 
Table  15  describes  the  four  scenarios  for  the  Elizabeth  River  study. 


Table  15.  Simulation  matrix  for  Elizabeth  River  using  WASP  7. 


13-box  Elizabeth 
River 

Hydrodynamic 

Dispersion 

Settling/Resuspension 

Velocities 

Mineralization 

Case  3a 

Yes 

No 

No 

Case  3b 

Yes 

No 

Yes 

Case  4a 

Yes 

Yes 

No 

Case  4b 

Yes 

Yes 

Yes 

Case  3a  simulates  the  transport  of  TNT  in  Elizabeth  River  with  the  8-kg/year  TNT  released  in 
Box  #2,  #6,  and  #12,  respectively.  The  goal  of  evaluating  the  behavior  at  these  three  different  release 
locations  is  to  understand  the  transport/dispersion  patterns  for  TNT  in  the  river. 

Case  3b  is  similar  to  Case  3a,  but  with  mineralization. 

Case  4a  is  similar  to  Case  2a  and  simulates  both  transport  in  water  column  and  exchange  of 
particulate  TNT  across  the  water/sediment  interface  for  Elizabeth  River. 

Case  4b  is  the  same  as  Case  2a,  but  with  mineralization. 

3.7  MODELING  PORE-WATER  ADVECTION  IN  SEDIMENT  FLUME  BOX 

As  part  of  the  empirical  study  for  evaluating  MC  release  in  sediment,  experiments  were 
conducted  to  measure  pore-water  advection  flow  and  release  rates  of  dye  (as  MC  surrogate)  at  the 
University  of  Wisconsin-Milwaukee.  Results  of  this  empirical  study  were  then  compared  with  the 
simulations  of  pore- water  advection.  To  simulate  the  pore-water  in  the  sediment  flume  box,  a  simple 
2-D  model  based  on  Darcy’s  Law  was  developed.  Figure  33  shows  the  configuration  of  the  model 
grid  and  parameterization.  Three  head  differences  between  the  left  and  right  boundaries  were  used:  1, 
2,  and  3  inches. 
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•  •  • 
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Figure  33.  2-D  model  for  pressure  induced  pore-water  seepage  flow  in  sand  flume  box. 

Following  Darcy’s  Law, 


Q  =  -kA—,  (32) 

as 

where  Q  is  the  Darey  veloeity  (L/sec);  k  is  the  conductivity  (L/sec),  which  is  a  function  of  the 
media  (e.g.,  silt,  sand,  etc);  A  is  the  cross  section  of  the  flow,  which  is  perpendicular  to  the  flow 
direction  5,  and  dh/ds  is  the  pressure  gradient  in  the  direction.  Eq.  (32)  is  a  well-known 
fundamental  governing  equation  for  groundwater  flow  and  transport,  which  is  a  basis  for  nearly 
all  groundwater  models  (MODFLOW,  SUTRA,  etc.)  (Bear,  1979  and  Freeze  and  Cherry,  1979). 
While  these  complex  models  can  accommodate  complex  systems,  the  fundamental  principal 
depicted  by  Eq.  (32)  is  universal. 

Similarly,  we  have  developed  a  simple  model,  based  on  Eq.  (32),  to  simulate  pore-water  advection 
for  the  same  scenarios  that  were  evaluated  by  the  empirical  study.  For  those  experiments,  it  is 
assumed  that  the  pore-water  flow  was  homogeneous,  isotropic  (no  directionality  for  conductivity), 
and  saturated  flows,  for  which  the  continuity  equation  is 

V.(pe)=0,  (33) 

where  the  gradient  vector  operator  (V)  and  the  flow  velocity  ( V)  can  be  expressed  in  the  (x,  y ,  z) 
coordinate  system  as 


and 


„  d  ~  d  ~  d  -r 

V  =  —  i  ^ - /  H - k 

dx  dy  dz 


V  =  VJ  +VJ  +V^k  . 


(34) 


(35) 
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Substituting  Eqs.  (32)  into  (33)  (Q  =  V),  the  equation  governing  the  homogeneous,  isotropic 
saturated  flow  for  the  empirical  sand  flume  box  is 


d^h  d^h 
dx^  ^  dy^ 


=  0. 


(36) 


Assuming  a  rectangular  grid  is  used  with  the  size  dx  =  dy=  3  inches,  60  cells,  each  with  the  I- 
index,  correspond  to  the  x-direction,  and  12  cells,  each  with  the  J-index,  correspond  to  the  y- 
direction.  The  model  calculates  the  head  pressure,  h(I,J),  induced  by  the  pressure  head  difference 
imposed  at  the  two  boundaries  (Figure  34).  The  indices,  (I,J),  represent  the  coordinates,  (x,  y)  with 
equal  grid  sizes,  dx  =  dy  =  3,  and  x  =  I*dx,  y  =  J*dy,  respectively.  In  order  to  solve  Eq.  (36)  at  all  the 
nodal  points  (I,J),  the  head  at  the  node  (I,J)  can  be  categorized  into  the  following  three  types;  interior 
nodes,  boundary  nodes,  and  comer  nodes  (Figure  34). 


J 


- 

- ir 

- 

I 


J 


Interior  node  Boundary  node  Corner  node 

Figure  34.  Interior,  boundary,  and  corner  nodes  of  the  numerical  mesh  for  solution  of  Eq.  (36). 

Boundary  and  comer  nodes  in  the  reversed  direction  can  be  developed  following  the  same  mles  as 
shown  in  the  following  equations. 

To  set  up  the  model  following  Eq.  (36),  the  equations  are  applied  for  the  nodes  in  the  interior,  at 
the  boundaries,  and  at  the  comers  of  the  domain; 

Interior  nodes;  h(I,  -  b  «^)  +  h{I  +  \,J)  +  h{I,J  -\)  +  h{I,  J  + 1)) 

Comer  nodes;  h(I,J)  =  ^{h(I,J  +  l)  +  h(I  + 1,  J)) 

The  model  calculates  the  pressure  head,  h(I,J),  at  each  of  the  nodes.  Figure  35  (top  panel)  shows 
the  pressure  head  and  the  corresponding  pressure  gradient  contours  for  the  1-inch  head  difference 
(Ah)  between  the  boundaries.  The  other  two  panels  in  the  figure  show  the  corresponding  variables  for 
the  other  two  cases,  2-inch  and  3-inch  head  difference,  respectively.  Darcy’s  velocity  can  be  readily 
obtained  as  the  product  of  the  head  gradients  and  the  conductivity  of  the  sand.  Generally,  pore-water 
advection  is  proportional  to  the  net  head  difference  between  the  two  boundaries.  Flows  (arrows)  are 
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perpendicular  to  the  head  contours  (solid  line  with  numbers)  flowing  from  high  head  (left)  to  low 
head  (right).  In  the  middle  section  of  the  box,  pore- water  flows  are  relatively  uniform  in  the  water 
column. 


Sand  Flume  Box  Width  (15  Inches) 

Figure  35.  Simulated  pressure  gradients  over  the  three  scenarios  (Ah  =  1,2,  and  3  inches,  for  the 
top,  middle,  and  bottom  figures,  respectively). 
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4.  RESULTS  AND  DISCUSSION 


4.1  MODEL  VALIDATION 

4.1.2  Validation  Case  I:  Dependence  on  Current  Speed,  (U),  and  Radial  Direction,  (r) 

In  the  first  validation  case,  the  radius  of  the  breach  hole  isb  =  5  mm.  The  inside  cavity  is  a  semi¬ 
sphere  with  radius  i?  =  0.1  m.  Upstream  current  speed  varies  from  1  to  0.5  m/s.  The  Reynolds 
number  based  on  the  hole  size  and  the  mean  current  speed  varies  from  5  to  2,500. 

The  flow  configuration  results  in  a  complicated  circulation  pattern  inside  the  shell  as  seen  in 
Figure  36.  Circulation  of  fluid  causes  advection  of  MC  solution,  and  the  resulting  effects  can  be 
modeled  as  a  diffusion  process  with  a  circulation  enhanced  diffusivity  Da,  and  solution  of  Eqs.  (10a) 
and  (10b)  can  still  apply. 


Figure  36.  Current-induced  inner  circulation  in  the  casing,  b  =  5mm  and  U  =  6  cm/s. 

The  circulation  enhanced  diffusivity.  Da,  has  a  dimension  of  [Velocity  x  Length],  which  in 
general,  should  be  a  function  of  r;  Da  varies  with  distance  r  inside  the  cavity.  The  velocity  scale  is 
modeled  as 


w^(r) 


I  dA 

•'=r 


(37) 


where  u,  v,  and  w  are  the  disturbed  velocity  components  in  the  x-y-z  coordinates  inside  the  cavity  at 
the  radius  r,  over  which  the  surface  integral  is  performed. 
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The  value  is  the  averaged  velocity  on  a  spherical  surface  with  radius  r.  Da  is  modeled  as 


Da  =Pu^r, 


(38) 


where  p  is  a  model  parameter. 

Based  on  FLUENT  simulation  results,  Da/P  as  a  function  of  r  is  calculated  for  different  current 
speed,  as  shown  in  Figure  37.  Generally,  ua  decreases  with  r,  but  Ua  r  is  only  weakly  dependent  on  r, 
except  near  the  breach  hole  and  the  solid-liquid  interface. 


Figure  37.  Distribution  of  UAr/Ub  as  a  function  of  r  for  different  current  speed. 

As  a  first-order  approximation,  we  assume  Da  is  homogenous  inside  the  breach  hole  and  we  take 
the  average  of  Da  /p  calculated  at  different  r  (except  for  r  =  h  and  r  =  R)  as  the  effective  advection 
induced  diffusivity,  as  defined  in  Eq.  (38),  and  plot  it  with  respect  to  current  speed  U  times  the  hole 
radius  b,  as  shown  in  Figure  38. 

Again,  as  a  first-order  approximation,  we  assume  there  is  a  linear  relationship  between  Da/P  and 
Ub.  By  linearly  fitting  with  the  simulation  result  (Figure  39),  we  can  propose  a  model  for  Da'. 


^  =  0.0l2Ub. 

A 


(39) 


It  can  be  shown  that  Da  is  much  greater  than  the  molecular  diffusivity  Dm  except  for  very  small 
current  speed.  According  to  this  argument,  Eq.  (1 1)  for  the  MC  flux  rate  is 


F 


iTixO.OU^UbC^ 
2x0.0l2/3Ub  T 
aUb^  b 


0u[7dj^)Cs, 


(40) 
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where  6  =  ,  i.e.,  the  flux  is  linearly  proportional  to  the  current  speed,  the  saturated 

0.024/?  +  a 

MC  concentration,  and  the  radius  of  the  breach  hole. 


Figure  38.  Modeled  relationship  between  the  effective  diffusivity  Da,  and  the  current  speed  U. 


The  above  model  is  validated  with  the  simulation  results,  as  shown  in  Figure  39.  The  linear 
relation  observed  above  fits  well  with  simulation  results  and  the  model  coefficient  6  is  ~0.023.  Note 
that  the  linear  growth  of  F  with  U  is  only  valid  for  large  Reynolds  number  values.  As  shown  in  the 
inset  plot  for  small  U  in  Figure  39,  the  effect  of  molecular  diffusion  is  on  the  same  order  of 
magnitude  as  advection,  thus  the  linear  dependency  is  not  valid  in  this  regime. 

With  the  fitted  parameter  0,  we  can  substitute  the  modeled  Da  into  the  equation  of  concentration 
profile  Eq.  (12)  and  compare  it  with  simulation  results.  Figure  40  shows,  as  a  function  of  r,  both  the 
analytical  solution  and  simulation  results  for  different  current  speeds.  There  are  generally  good 
agreements,  especially  for  smaller  current  speed.  However,  for  D  >  100  mm/s  (+  symbols), 
differences  between  simulation  and  analytical  solutions  become  apparent.  This  is  due,  in  part,  to  the 
circulation-induced  diffusivity  becoming  increasingly  inhomogeneous,  which  makes  our  analytical 
solution  invalid  for  this  current  regime. 

4.1.3  Validation  Case  II:  Dependence  on  Hole  Size,  (b) 

As  in  Validation  Case  I,  the  cavity  inside  the  shell  casing  is  still  a  semi-sphere  with  radius 
i?  =  0.1  m  and  the  current  speed  is  a  constant  with  D  =  0.1  m/s.  The  radius  of  breach  hole  varies  from 
Z?  =  5  to  50  mm. 

If  we  first  assume  our  previous  analysis  is  still  valid  for  this  case,  then  we  would  have 

—  «  0.023D(;rZ)")=  0.007h^  (41) 

where  F  would  be  a  quadratic  function  of  the  hole  radius  b  (or  alternatively  a  linear  function  of 
breach  area  A).  Using  these  two  validation  cases,  the  MC  release  rate  is  approximated  as 


F  ^OmZUACs,  (42) 

similar  to  the  analysis  of  Eq.  (12)  in  Figure  39. 
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The  above  model  is  compared  with  the  simulation  results  in  Figure  41.  It  agrees  well  with 
numerical  experiment  results  for  b  <  25  mm.  However,  it  does  over-estimate  the  flux  for  larger  holes. 
This  is  not  a  surprise  since  when  b  grows  larger,  the  inner  circulation  becomes  more  complex  and  the 
assumption  of  a  small  hole  size  starts  to  be  violated  for  the  analytical  release  rate  function. 


Figure  39.  Variation  of  MC  flux  rate  with  respect  to  current  speed. 
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Figure  40.  Radial  distribution  of  concentration  in  the  casing.  Simulation  results  compared  with 
analytical  solution  (Eq.  (9)). 
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Figure  41 .  Variation  of  MC  flux  rate  F,  as  a  function  of  hole  radius  b. 
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4.2  EFFECTS  OF  SHELL  THICKNESS 

In  our  previous  analysis,  the  munition  wall  or  shell  thickness  (hereafter  denoted  as  h)  has  been 
ignored.  Shell  thickness  is  generally  smaller  compared  to  the  inner  cavity  radius  R,  thus  its  effects  on 
the  overall  flux  rate  can  be  considered  secondary,  i.e.,  it  may  result  in  a  slightly  different  coefficient 
P  in  the  boundary  condition  of  our  analytical  model  {i.e.,  Eq.  (38)),  but  the  general  equation  of  inner 
scalar  transport  should  not  be  affected  by  nominal  shell  thickness. 

Due  to  the  limitation  of  grid  size,  details  of  circulation  patterns  smaller  than  h  are  not  easily 
resolved  by  simulation.  We  will,  however,  present  the  empirical  relation  between  the  release  rate  and 
the  shell  thickness.  Three  breach  radius  are  selected  for  simulation,  b  =  2.5,  5,  and  10  mm, 
respectively.  The  shell  thickness  varies  from  4  to  24  mm  and  the  current  velocity  is  allowed  to 
change  from  4  to  45  cm/s.  Figure  42  shows  the  release  rate  function  F  v^'.  the  ambient  current  U  for 
various  shell  thicknesses.  In  most  cases,  a  linear  increase  of  release  rate  (flux)  with  the  mean  current 
U  can  be  observed.  This  agrees  with  the  analytical  solution  Eq.  (10a).  However,  the  proportionality 
{6X)  differs  from  case  to  case.  For  comparison,  the  empirical  solution  for  validation  case  I,  Eq.  (40), 
is  also  shown  in  this  figure. 

The  release  rate  function  (Eq.  (42))  is  shown  normalized  and  plotted  v^'.  shell  thickness  in  Figure 
43.  The  normalized  release  rate  fluctuates  when  h  and  b  are  of  the  same  order  of  magnitude.  When  h 
»  b,  the  normalized  F  seems  to  be  independent  of  h  and  becomes  a  constant.  This  constant  varies 
between  0.015  and  -0.035,  compared  to  0.023  demonstrated  in  the  previous  validation  case. 
Therefore,  future  efforts  will  select  an  appropriate  value  of  this  constant  based  on  the  h:b  ratio  that 
lies  within  the  range  of  0.015  to  -0.035. 

4.3  RELEASE  TIME 

The  general  analytical  release  rate  function  (Eq.  (10a)),  which  has  now  been  validated  by  the 
FLUENT  model,  can  be  used  to  predict  the  release  rates  of  MC,  under  various  hydrodynamic  and 
shell  integrity  conditions.  For  a  single  shell  with  known  MC  mass,  the  predicted  release  rates  can  be 
further  used  to  estimate  time  to  deplete  all  MC  in  the  shell  under  these  conditions.  Here  we  use  Eq. 
(10a)  for  the  calculations  of  release  time  for  TNT,  RDX,  and  HMX  under  various  hydrodynamic  and 
shell  integrity  conditions. 

In  the  analytical  release  rate  function,  the  MC-dependent  variables  include  the  saturation 
concentrations  (Cs)  and  the  dissolution  speed  (p),  obtained  by  dividing  the  dissolution  rate  by  the 
saturation  concentration.  These  MC-dependent  chemical  data  are  listed  in  Table  16. 

Brannon  et  al.  (2005)  measured  dissolution  rate  constants  of  TNT,  RDX,  and  HMX  in  saline  water 
and  compared  them  with  the  rate  constants  in  freshwater  (Table  16). 

“Dissolution  rates  of  RDX  and  HMX from  the  military  formulations,  Composition  B  and 
Octol,  did  not  differ  significantly  under  freshwater  and  saline  conditions  (see  Figure  44).  The 
rate  of  dissolution  of  TNT from  the  solid  in  saline  water  was  slower  than  in  freshwater. 
However,  freshwater  dissolution  rate  values  for  pure  TNT  would  conservatively  estimate 
TNT  saline  water  dissolution  rates.  The  rate  of  dissolution  of  TNT  at  salinities  of  10  and 
30  ppt  was  the  same  as  the  rate  at  20  ppt.  Dissolution  rates  for  the  explosives  in  both  fresh 
and  saline  water  decreased  in  the  order  TNT>RDX>HMX,  which  is  consistent  with  trends 
observed  previously  in  fresh  water...  ”  Brannon  et  al.  (2005^ 
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b  =  2.5  mm 


U(m/s) 
b  =10.0  mm 


Figure  42.  Release  rate  vs.  current  speed  for  selected  shell  thickness,  (h),  and  breach  hole  radii,  (b). 
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Figure  43.  Effects  of  shell  thickness  on  the  normalized  release  rate. 


Table  16.  Saturation  concentrations  (Cs),  dissolution  rates,  and  dissolution  speed  (p)  for  TNT,  RDX, 
and  HMX  in  freshwater. 


MCs 

Saturation 
Concentration  (Cs) 
(mg/L) 

Dissoiution  Rate 
(jjg/cm^/s) 

Dissoiution  Speed  (p)  (cm/s) 
(Dissoiution  Rate/Saturation 
Concentration) 

TNT 

88.5  (a) 

0.20(d) 

0.00226 

RDX 

38.4  (b) 

0.05(d) 

0.00130 

HMX 

6.6  (c) 

0.15(d) 

0.02272 

(a)  Ro  et  al.  (1996) 

(b)  Spanggord  et  al.  (1983) 

(c)  McLellan,  Hartley,  and  Brower  (1988) 

(d)  Lynch,  Brannon,  and  Delfino,  2002a  and  2002b. 

Therefore,  the  dissolution  rates  for  freshwater  (Table  16)  are  used  as  a  conservative  measure  for 
the  analysis  of  release  rate  of  MCs  from  breached  shells  in  the  following  sections. 

For  analysis  of  the  following  scenario,  we  use  a  typical  round  with  the  following  shell  dimensions: 
length  =  24  inches,  diameter  =  8  inches,  and  MC  weight  =  8  kg.  Shell  thickness  is  not  included  for 
problem  simplification  purposes. 

For  different  realistic  combinations  of  the  five  variables  in  Eq.  (10a),  release  time  can  vary  by 
several  orders  of  magnitude.  For  example,  under  an  ambient  current  of  45  cm/s  (fast  flow),  release 
time  ranges  from  300  years  for  a  hole  with  radius  of  0.05  cm  to  1  year  for  a  hole  with  radius  of  0.95 
cm,  illustrating  the  importance  of  hole  size.  The  sensitivity  to  ambient  current  is  also  apparent  in 
Figure  45,  where  the  release  time  varies  by  one  order  of  magnitude  between  slow  flow  (5  cm/s)  and 
fast  flow  (45  cm/s).  For  example,  for  a  hole  with  a  radius  of  0.5  cm,  it  takes  about  28  years  to  deplete 
the  MC  inside  the  shell  under  ambient  current  of  5  cm/s.  The  release  time  is  much  shorter,  less  than 
8  years,  if  the  ambient  current  increases  to  45  cm/s. 
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Figure  46  compares  the  release  time  for  TNT,  RDX,  and  HMX  under  the  same  conditions.  TNT 
release  time  is  fast,  followed  by  RDX,  and  HMX.  Although  HMX  has  the  highest  dissolution  speed, 
its  low  saturation  concentration  (Table  16)  and  interactions  with  other  parameters,  including  current 
speed,  hole  size,  and  cavity  size,  make  the  mass  release  from  the  shell  the  slowest  of  the  three  MCs 
for  this  scenario.  For  the  scenario  in  Figure  46,  saturation  concentrations  dictate  the  release  time  (and 
rate)  for  the  three  MCs. 

In  Figure  47,  time  for  complete  release  of  MC  from  a  hypothetical  shell  depends  on  the  five 
parameters,  as  discussed  previously  with  Eq.  (10a).  Overall,  the  release  function  is  controlled  by  the 
dissolution  rate  for  the  regime,  pR^  «  aUb^.  In  this  region,  the  release  rate  (and  the  time  for 
complete  release)  depends  primarily  on  the  dissolution  rate  (speed),  but  is  insensitive  to  the  current, 
which  is  assumed  to  be  large  enough  to  disperse  released  solution  faster  than  the  dissolution  rate.  For 
the  other  extreme,  the  current  controlled  flow  regime,  pR^  »  aUb^,  the  release  function  (and  the 
time  for  complete  release)  is  less  dependent  on  the  dissolution  rate,  and  depends  primarily  on  the 
current  speed.  The  dissolution  speeds  are  0.0016  cm/s  for  TNT,  0.0008  cni/s  for  RDX,  and  0.019 
cm/s  for  HMX,  respectively,  which  the  three  up-point  arrows  indicate  as  release  processes  for  the 
three  MCs  are  between  the  two  limiting  cases. 


Figure  44.  Measured  dissolution  rate  constants  for  TNT,  RDX,  and  HMX  in  freshwater  and  saline 
water  (Brannon  et  al.,  2005). 

In  a  natural  coastal  setting,  tidal  cycles  dictate  the  flow  behavior.  Release  rates  during  tidal  cycles 
with  tidal  current  speed  ranging  from  5to  45  cm/s  are  shown  in  Figure  48  and  Figure  49.  For  small 
hole  size  (b  =  0.1  cm.  Figure  48),  release  rates  for  TNT,  RDX,  and  HMX  vary  with  tidal  currents 
over  the  entire  current  speed  range.  For  this  scenario,  the  current-based  Reynolds  Number  (Ub/D)  is 
much  smaller  than  the  dissolution-based  Reynolds  Number  (pR/D)  and  release  rates  are  limited  by 
current  speed.  For  larger  holes  {i.e.,  b  =  2.5  cm.  Figure  49),  release  rate  of  TNT  is  limited  primarily 
by  current  speed  when  the  current  is  low  (<  15  cm/s).  When  the  current  becomes  strong  (speed  ~15 
to  45  cm/s),  release  rate  of  TNT  is  only  weakly  dependent  on  current  speed;  it  is  limited  primarily  by 
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dissolution  rate.  For  RDX,  Figure  48  and  Figure  49  show  that  release  rate  is  primarily  limited  by 
current  speed  for  small  hole  size  {i.e.,  b  =  0.1  cm),  but  it  is  primarily  limited  by  dissolution  rate  for 
larger  hole  size  {i.e.,  b  =  2.5  cm)  over  the  entire  range  of  current  speed.  The  release  rate  for  HMX  is 
limited  by  current  speed  for  both  small  and  large  hole  sizes  over  the  entire  range  of  current  speed. 

From  these  analyses,  it  has  been  shown  that  release  rates  for  TNT,  RDX,  and  HMX  are  governed 
by  five  variables,  including  current  speed,  hole  size,  cavity  volume,  dissolution  rates,  and  hydro- 
dynamic  diffusivity  coefficient.  Release  rates  of  MC  in  a  dynamic  estuarine  environment  are  limited 
by  two  overall  conditions;  current-controlled  release  and  dissolution-controlled  release.  Ultimately, 
the  release  rate  function  can  be  linked  with  a  hydrodynamic  and  fate/transport  model  to  simulate  the 
dynamic  release  and  subsequent  fate/transport  of  MCs  in  the  estuarine  receiving  water. 

TNT  Release  Time  (L=M'',  8  Kg) 
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Figure  45.  Days  to  complete  release  of  TNT  inside  the  shell  with  dimension  of  L  =  24  inches,  D  =  8 
inches,  and  MC  mass  =  8  kg.  Release  time  are  plotted  as  functions  of  ambient  velocity  (U),  and  hole 
radius  (b). 


TNT-RDX-HMX 
(Shell:  L=24",  D=8",  8kg) 


Figure  46.  Release  time  for  TNT,  RDX,  and  HMX  under  the  same  conditions  (U  =  15  cm/s). 
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Figure  47.  Days  for  complete  release  (based  on  Eq.  (10a))  for  TNT,  RDX,  and  HMX  for  a  release 
case:  U  =  10  cm/s,  R  =  10  cm.  The  three  up-pointing  arrows  represent  the  dissolution  speeds  for  the 
three  compounds. 


Figure  48.  Release  rates  and  current  speeds  over  two  tidal  cycles  (~24  hours)  for  TNT,  RDX,  and 
HMX  (b  =  0.1  cm,  R=  10  cm). 
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Figure  49.  Release  rates  and  current  speeds  over  two  tidal  cycles  (~24  hours)  for  TNT,  RDX,  and 
HMX  (b  =  2.5  cm,  R  =  10  cm). 

4.4  MODELING  FATE  AND  TRANSPORT  OF  TNT,  RDX,  AND  HMX 
4.4.1  TRIM2D  Model  Results 

Source  release  rates  are  a  function  of  current  speed  as  determined  by  tidal  fluctuations,  as  illus¬ 
trated  in  Figure  50.  Current  speeds  in  the  inner  (San  Diego)  bay  area  are  on  the  order  of  8  cm/s  and 
the  corresponding  release  rates  are  about  300  pg/s.  Generally,  the  release  rate  is  controlled  by  current 
speed,  below  current  speeds  of  12  cm/s.  When  current  speeds  exceed  12  cm/s,  the  release  rate  is 
limited  by  the  dissolution  rate  and  becomes  insensitive  to  a  high  current  speed.  Release  rates  under 
these  extreme  conditions  (high  current  speed  and  ignoring  degradation  processes)  were  previously 
analyzed  analytically,  and  they  are  now  demonstrated  here  under  field-based  conditions. 

Scenario  A:  Steady-state  TNT  concentrations  in  water  from  a  single  shell  source  release  simulated 
under  various  tidal  stages  without  any  degradation  are  shown  in  Figure  5 1  and  Figure  52.  Generally, 
steady-state  concentrations  decrease  from  the  source  to  the  mouth  of  the  bay,  resulting  from 
hydrodynamic  transport  and  diffusivity.  Under  steady-state  conditions,  transport  of  TNT  released 
from  the  breached  shell  gradually  achieves  equilibrium  with  the  hydrodynamic  flushing  of  the  bay 
water,  which  is  driven  by  ocean  tides  mixed  with  diurnal  and  semi-diurnal  tidal  cycles.  The  simulated 
steady-state  TNT  concentrations  range  from  ~30  to  50  ng/L  in  the  inner  bay  to  10  to  30  ng/L  in  the 
mid-bay  areas  and  to  about  1  to  2  ng/L  near  the  bay  mouth. 

Results  shown  in  Figure  52  and  Figure  53  are  independent  of  dissolution  processes  and  illustrate 
only  the  transport  behavior  of  TNT;  it  is  assumed  that  already  dissolved  TNT  (i.e.,  to  simulate 
situation  where  dissolution  rate  is  not  limiting)  is  released  into  the  water  column.  No  net  loss  of  TNT 
from  the  water  column  is  allowed  to  exist  from  any  other  process,  except  the  hydrodynamic 
dispersion  and  tidal  flushing  through  the  mouth  of  the  bay.  Hydrodynamic  dispersion  redistributes 
TNT  concentrations  and  tidal  flushing  through  the  mouth  of  the  bay  provides  a  net  mass  loss  of  TNT 
out  of  the  model  domain.  When  a  steady-state  condition  is  reached,  mass  loss  rate  of  TNT  by 
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flushing  through  the  mouth  of  the  bay  should  be  balanced  by  the  net  release  rate  of  TNT  from  the 
breached  shell.  Therefore,  results  shown  in  Figure  52  and  Figure  53  can  be  considered  conservative 
estimates  of  TNT  exposure  from  single-shell  source  releases. 

Release  Rates  (Inner  Bay) 
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- Release  Rate  (ug/s)  ^—Current  Speed  (cm/s) 


Figure  50.  Current  speed  and  estimated  release  rate  at  the  inner  San  Diego  Bay  location. 
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Figure  51 .  Simulated  TNT  concentration  at  the  source  location  (inner  San  Diego  Bay)  and  total  TNT 
load  during  1-year  period. 
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Figure  52.  Steady-state  TNT  concentration  during  ebb  tide  (simulated  release  from  inner  San  Diego 
Bay)  (x-y  coordinates  are  distances  with  I,  J  indices  multiplied  by  100  meters,  the  same  convention 
hereafter). 


In  a  separate  study,  the  mineralization  of  TNT,  RDX,  and  HMX  in  marine  water  and  sediments 
was  reported  in  Montgomery  et  al,  2008. 

These  mineralization  rates  represent  neither  a  zero  or  first-order  process.  For  modeling  purposes, 
we  used  these  reported  mineralization  measurements  to  estimate  and  incorporate  pseudo-first-order 
TNT  degradation  kinetics  into  the  models.  As  shown  in  Table  17,  1.5%  of  the  TNT  is  mineralized  in 
marine  sediment  and  seawater  by  the  end  of  the  first  day.  This  daily  mineralization  rate  is  equivalent 
to  a  first-order  mineralization  rate  of  0.0151  day  ^  Simulations  with  TRIM2D  using  this  mineraliza¬ 
tion  rate  were  conducted  and  the  resulting  simulated  steady-state  concentrations  are  shown  in  Figure 
54  and  can  be  compared  to  the  results  without  mineralization  as  shown  in  Figure  52. 

Scenario  B.  For  this  scenario,  current  speeds  near  the  Naval  Station  San  Diego  are  on  the  order  of 
4  cm/s  and  the  corresponding  release  rates  are  approximately  150  ug/s.  Release  rate  during  the  entire 
1-year  period  is  current-controlled,  due  to  the  relative  low  current  speed  and  tidal  fluctuations,  as 
shown  in  Figure  55. 

Simulated  steady-state  TNT  concentrations  are  shown  in  Figure  56  through  Figure  58.  Generally, 
as  observed  previously  in  Scenario  #1,  steady-state  concentrations  decrease  from  the  source, 
resulting  from  hydrodynamic  transport  and  dispersion  within  the  bay.  Simulated  steady-state  TNT 
concentration  ranges  were  ~8  ng/L  in  the  inner  bay,  ~4  ng/L  in  the  mid-bay  areas,  and  ~1  ng/L  near 
the  bay  mouth.  The  highest  TNT  concentrations  (~15  ng/L)  occur  near  the  source  region,  at  Naval 
Station  San  Diego.  Due  to  the  weaker  current  speed  at  this  location  in  the  bay,  annual  release  of  TNT 
from  the  breached  shell  at  Naval  Station  San  Diego  is  only  about  5  kg/year,  in  comparison  with  the 
larger  release  rate  of  about  8.4  kg/year  in  the  southern  inner  bay  (Scenario  #1)  (Figure  51).  Simulated 
steady-state  concentrations  taking  into  account  the  previously  described  first-order  mineralization 
rate  are  shown  in  Figure  59. 
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In  a  separate  study,  the  mineralization  of  TNT,  RDX,  and  HMX  in  marine  water  and  sediments 
was  reported  in  Montgomery  et  al,  2008. 

These  mineralization  rates  represent  neither  a  zero  or  first-order  process.  For  modeling  purposes, 
we  used  these  reported  mineralization  measurements  to  estimate  and  incorporate  pseudo-first-order 
TNT  degradation  kinetics  into  the  models.  As  shown  in  Table  17,  1.5%  of  the  TNT  is  mineralized  in 
marine  sediment  and  seawater  by  the  end  of  the  first  day.  This  daily  mineralization  rate  is  equivalent 
to  a  first-order  mineralization  rate  of  0.0151  day  '.  Simulations  with  TRIM2D  using  this  mineraliza¬ 
tion  rate  were  conducted  and  the  resulting  simulated  steady-state  concentrations  are  shown  in  Figure 
54  and  can  be  compared  to  the  results  without  mineralization  as  shown  in  Figure  52. 

Scenario  B.  For  this  scenario,  current  speeds  near  the  Naval  Station  San  Diego  are  on  the  order  of 
4  cm/s  and  the  corresponding  release  rates  are  approximately  150  ug/s.  Release  rate  during  the  entire 
1-year  period  is  current-controlled,  due  to  the  relative  low  current  speed  and  tidal  fluctuations,  as 
shown  in  Figure  55. 

Simulated  steady-state  TNT  concentrations  are  shown  in  Figure  56  through  Figure  58.  Generally, 
as  observed  previously  in  Scenario  #1,  steady-state  concentrations  decrease  from  the  source, 
resulting  from  hydrodynamic  transport  and  dispersion  within  the  bay.  Simulated  steady-state  TNT 
concentration  ranges  were  ~8  ng/L  in  the  inner  bay,  ~4  ng/L  in  the  mid-bay  areas,  and  ~1  ng/L  near 
the  bay  mouth.  The  highest  TNT  concentrations  (~15  ng/L)  occur  near  the  source  region,  at  Naval 
Station  San  Diego.  Due  to  the  weaker  current  speed  at  this  location  in  the  bay,  annual  release  of  TNT 
from  the  breached  shell  at  Naval  Station  San  Diego  is  only  about  5  kg/year,  in  comparison  with  the 
larger  release  rate  of  about  8.4  kg/year  in  the  southern  inner  bay  (Scenario  #1)  (Figure  51).  Simulated 
steady-state  concentrations  taking  into  account  the  previously  described  first-order  mineralization 
rate  are  shown  in  Figure  59. 
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Figure  53.  Steady-state  TNT  concentration  during  flooding  tide  (simulated  release  from  inner 
San  Diego  Bay). 
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Figure  54.  Steady-state  TNT  concentrations  with  first-order  mineralization  during  ebb  tide  (simulated 
release  from  inner  San  Diego  Bay). 
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Table  17.  Mineralization  rates  for  TNT  (data  provided  by  Dr.  Mike  Montgomery,  NRL). 


Microbial  Assemblage  Type 
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Figure  55.  Current  speed  and  estimated  release  rate  at  Naval  Station  San  Diego. 
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Figure  56.  Simulated  TNT  concentration  at  the  source  location  (Naval  Station  San  Diego)  and  total 
TNT  load  during  1-year  period. 
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Figure  57.  Steady-state  TNT  concentrations  during  ebb  tide  (simulated  release  from  Naval  Station 
San  Diego). 
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Figure  58.  Steady-state  TNT  concentrations  during  flooding  tide  (simulated  release  from  Naval 
Station  San  Diego). 
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Figure  59.  Steady-state  TNT  concentrations  with  mineralization  of  TNT  during  flooding  tide 
(simulated  release  from  Naval  Station  San  Diego). 
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For  these  simulations  with  TRIM2D,  as  with  Scenario  #1,  it  is  assumed  that  all  TNT  released  from 
the  breached  shells  is  in  the  form  of  dissolved  species  (aqueous  TNT  concentrations).  The  only  TNT 
mass  loss  term  included  is  a  result  of  flushing  dissolved  TNT  out  of  the  bay  by  currents  and 
hydrodynamic  dispersion.  No  loss  to  the  sediment  is  assumed.  Therefore,  this  approach  focuses  on  a 
maximum  dissolved  TNT  concentration,  and  fine-resolution  model  results  from  TRIM2D  can  be 
considered  upper-bound  estimates  for  steady-state  TNT  concentrations  in  the  water  column.  Thus, 
the  predicted  steady-state  TNT  concentrations  using  TRIM2D  throughout  the  bay  can  be  considered 
conservative  estimates. 

4.4.2  WASP7  Model  Results 

For  verification  for  the  San  Diego  Bay  model,  WASPV-simulated  results  are  compared  with  those 
previously  simulated  by  TR1M2D.  To  compare  with  results  from  the  TR1M2D  model,  which  are 
considered  conservative  predictions  for  water  column  concentration,  a  steady-state  TNT  loading  rate 
of  8.4  kg/year  (Case  #1  for  TRIM2D  modeling  scenario)  was  also  used  for  loading  into  WASP?  Box 
#1  (refer  to  Figure  30).  Model  simulations  were  run  until  model  results  reached  steady-state.  Figure 
60  shows  the  simulated  steady-state  TNT  concentrations  for  the  eight  boxes  of  the  water  column  of 
the  San  Diego  Bay  for  Cases  la  and  lb.  For  Case  la,  hydrodynamic  dispersion  is  the  only  driving 
force  for  the  TNT  transport.  Simulated  TNT  concentrations  in  the  water  column  range  from  ~45  ng/L 
in  the  inner  bay  (Box  #1),  where  the  TNT  source  is  released,  to  ~5  ng/L  in  the  middle  bay  and  <  1 
ng/L  near  the  mouth.  TRIM2D  results  are  volume-averaged  for  each  of  the  eight  boxes  and  compared 
with  WASP  results.  As  shown  in  Figure  60,  the  WASP-simulated  TNT  concentrations  are  very  close 
to  the  TRIM-simulated  TNT  concentrations,  even  though  these  two  models  use  different  solution 
schemes  and  grid  resolutions.  The  results  for  Case  lb  also  compare  well  with  the  results  from 
TRIM2D  results  (Figure  60)  for  the  corresponding  scenarios. 
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Figure  60.  Wasp7-simulated  TNT  concentrations  in  water  column  for  Cases  1a  and  1b. 
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Unlike  TRIM2D,  WASP7  can  simulate  exchange  of  solids  and  particulate  contaminants  across  the 
water/sediment  interface  by  partitioning  and  settling/resuspension.  Two  partitioning  coefficients,  0.9 
and  5.0  L/Kg,  represent  the  lower  and  upper  bounds  from  the  literature  values  (Case  2a).  As 
discussed  in  the  literature,  partitioning  coefficients  for  TNT  (in  soil)  are  low  and  hence,  we  would 
expect  that  only  a  small  portion  of  TNT  mass  would  be  trapped  in  sediment.  As  shown  in  Figure  61, 
with  the  5.5  L/Kg  partitioning  coefficient,  TNT  concentrations  in  sediment  range  from  0.1  ng/L  for 
Box  #1  to  0.0001  ng/L  for  Box  #8.  With  a  smaller  partitioning  coefficient,  simulated  TNT 
concentrations  in  sediment  are  even  lower,  ranging  from  0.02  ng/L  for  Box  #1  to  0.00001  ng/L  for 
Box  #8.  While  these  concentrations  are  low,  note  that  these  are  simulated  concentrations  from  TNT 
release  from  a  single  breached  shell.  When  multiple  shells  exist  and  are  taken  into  account,  their 
accumulated  concentration  needs  further  study. 

Model  simulations  with  two  partitioning  coefficients,  i.e.,  5.54  and  0.88  L/kg,  with/without  a  first- 
order  mineralization  rate  constant  (0.0151  day'')  were  then  conducted.  Steady-state  TNT  concentra¬ 
tions  are  shown  in  Figure  61.  Generally,  TNT  concentrations  in  San  Diego  Bay,  which  is  similar  to 
most  estuaries,  are  primarily  governed  by  hydrodynamic  transport,  including  bay  flushing  and  hydro- 
dynamic  dispersion.  For  San  Diego  Bay,  hydrodynamic  transport  governs  the  primary  distributions 
of  TNT  released  from  a  single  breached  shell.  Steady-state  concentrations  result  from  an  equilibrium 
achieved  between  the  TNT  source  and  the  hydrodynamic  dispersion.  Generally,  TNT  concentrations 
in  sediment  are  proportional  to  a  partitioning  coefficient,  as  shown  in  Figure  61.  With  the  first-order 
mineralization,  steady-state  concentrations  in  sediment  decrease  by  20%  or  less. 
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Figure  61 .  WASP7-simulated  TNT  concentrations  in  sediment  for  Case  2a  (water  column  transport 
and  exchange  with  sediment)  and  2b  (with  mineralization)  for  different  partitioning  coefficients  for 
San  Diego  Bay. 
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For  the  Elizabeth  River  model  (Figure  29),  we  release  the  8  kg/year  TNT  at  three  different  boxes, 
#2,  #6,  and  #12  which  are  located  upstream,  middle  portion,  and  downstream  of  the  main  river 
stretch.  Unlike  San  Diego  Bay,  riverine  freshwater  in  Elizabeth  River  flows  from  upstream  all  the 
way  to  the  downstream  mouth  regions.  Therefore,  freshwater  flushing  exists  throughout  the  river. 
Figure  62  shows  that  release  in  Box  #2  generates  a  maximum  steady-state  concentration  of  ~42  ng/L. 
Release  in  Box  #6  produces  a  maximum  concentration  of  ~25  ng/L  at  Box  #6,  which  is  almost  equal 
to  the  concentration  resulting  at  Box  #6  upon  release  in  Box  #2.  The  same  phenomenon  is  observed 
for  release  in  Box  #12.  Release  in  Box  #12  produces  a  maximum  concentration  of  ~5  ng/L  at  Box 
#12,  which  is  close  to  the  concentration  resulting  at  Box  #12  upon  release  in  Box  #2.  These 
phenomena  are  consistent  with  processes  occurring  in  a  conservative  system.  Since  in  this  case,  there 
is  no  mineralization  and  the  transport  equation  satisfies  conservation  of  mass  law,  the  patterns  of  the 
predicted  concentrations  are  reasonable. 

With  the  mineralization  rate,  simulated  TNT  concentrations  in  the  water  column  show  reduced 
concentrations  in  the  water  column  as  shown  in  Figure  63.  On  average,  a  decrease  of -15%  in 
simulated  concentrations  is  observed. 

Simulated  TNT  concentrations  in  sediment  for  two  different  partitioning  coefficients,  0.88  and 
5.54  L/Kg  ,  both  with  and  without  mineralization  are  shown  in  Figure  64  (TNT  release  in  Box  #2) 
and  Figure  65  (TNT  release  in  Box  #6).  Simulated  TNT  concentrations  in  sediment  are  approximate¬ 
ly  proportional  to  the  partitioning  coefficient,  i.e.,  the  higher  partitioning  coefficients  are  associated 
with  higher  concentrations  in  sediment,  similar  to  the  results  for  San  Diego  Bay. 


Box  1  Box  2  Box  3  Box  4  Box  5  Box  6  Box  7  Box  8  Box  9  Box10  Box11  Box12  Box13 


Box  No. 


Figure  62.  Simulated  TNT  concentrations  in  water  column  for  Case  3a,  with  TNT  release  occurring 
in  Box  #2,  #6,  and  #12,  respectively. 
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Figure  63.  Simulated  TNT  concentrations  in  water  column  for  Case  3b,  with  TNT  release  occurring 
in  Box  #2,  #6,  and  #12,  respectively. 


Figure  64.  Simulated  TNT  concentrations  in  sediment  for  Case  4a  (with  partitioning  coefficients 
without  mineralization)  and  Case  4b  (with  mineralization),  with  TNT  release  occurring  in  Box  #2. 
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Figure  65.  Simulated  TNT  concentrations  in  sediment  for  Case  4a  (with  partitioning  coefficients  and 
without  mineralization)  and  Case  4b  (with  mineralization),  with  TNT  release  occurring  in  Box  #6. 

4.5  RESULTS  OF  TASKS  4.3.1A,  4.3.1B,  4.3.2,  AND  4.3.6 

Study  results  of  the  four  Phase  III  tasks  (Figure  3)  are  diseussed.  These  results  include  both  the 
numerical  modeling  study  and  empirical  study.  These  are  explained  as  follows. 

4.5.1  Effect  of  Shapes  of  Opening 

The  numerical  model,  FLUENT,  was  set  up  to  simulate  release  of  MC  from  holes  of  four  different 
shapes,  including  a  circle,  a  triangle,  and  two  rectangles  with  the  aspect  ratios  of  4:1  and  10:1, 
respectively.  Areas  for  these  four  different  opening  shapes  were  the  same,  which  is  equal  to  71/I6 
(cm^),  e.g.,  the  radius  of  the  circle  =  0.25  cm.  Simulated  release  rates  were  plotted  and  compared,  as 
shown  in  Figure  66. 

As  shown  in  Figure  66,  the  mass  release  rates  are  linearly  proportional  to  the  current  speed  (range 
of  -0.05-50  cm/s  (0.0005-0.5  m/s)),  with  two  relatively  constant  slopes  in  two  regions.  For  the 
region,  U  <  15  cm/s,  the  mass  release  rates  remain  almost  the  same  for  the  four  opening  shapes 
(Figure  67).  The  slope  is  about  0.0079  for  all  the  four  shapes  within  the  entire  current  speed  range 
(except  for  circle  at  0.15  m/s).  This  indicates  that  for  weak  flows,  shapes  of  the  hole  openings  do  not 
have  a  major  effect  on  the  MC  release  rate.  For  the  medium-high  flow  regimes,  namely,  U  >  0.15 
m/s,  the  mass  release  rates  are  linearly  dependent  on  the  current  speed,  but  with  different  slopes  for 
the  different  shapes  (Figure  68). 

Model  results  show  that  when  ambient  currents  are  weak  (less  than  0.15  m/s),  hydrodynamic 
mixing  of  MC  through  the  opening  is  also  weak  and  the  mass  release  rates  can  not  differentiate  the 
different  shapes  of  the  openings  and  differences  of  the  mass  release  rates  are  small. 
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Figure  66.  Four  different  shapes  of  opening  (above)  and  normalized  mass  release  rates  at  different 
current  speed  (lower). 

When  ambient  currents  become  stronger  (U  >  0.15  m/s),  hydrodynamic  mixing  of  MC  through  the 
opening  becomes  stronger  and  the  release  rates  are  influenced  by  the  shapes  of  the  openings.  Of  the 
four  different  shapes,  the  circle  produces  the  largest  release  rate  (with  a  slope  of  0.039),  followed  by 
the  rectangle  of  4: 1  aspect  ratio  (0.028),  and  then  by  the  triangle  (0.02)  and  the  rectangle  of  10: 1 
aspect  ratio  produces  the  smallest  mass  release  rate  (with  a  slope  of  0.01 1)  (Figure  68). 

4.5.2  Effect  of  Orientation  of  Opening  (rectangle  of  4:1  aspect  ratio) 

Breached  holes  are  formed  not  only  in  different  shapes,  but  also  at  different  orientation  angles 
relative  to  the  current  direction.  For  example,  breached  cracks  may  form  that  are  in  line  or  crossing 
with  the  ambient  current.  Effect  due  to  orientation  on  the  release  rate  was  studied  by  simulating 
release  rates  for  the  rectangle  (4:1)  with  the  orientation  angles  ranging  from  0°  (current  crosses  with 
the  long  axis),  to  45°  and  90°  (current  in  line  with  the  long  axis).  Figure  69  shows  that  a  45°-angle 
produces  the  maximum  release  rate  (F  =  0.041U-0.0033),  and  a  0°-angle  produces  the  minimum 
release  rate  (F  =  0.028U-0.004).  Such  behavior  of  different  orientation  angles  on  the  release  rate  may 
not  be  obvious  intuitively.  It  may  imply  that  local  dispersion  near  the  hole  at  different  orientation 
angles  might  be  more  complex  than  originally  thought.  The  model  results  show  that  differences  on 
the  release  rates  from  different  orientation  angles  are  in  the  range  of  +25%,  compared  to  the 
difference  of  twofold  on  the  release  rates  from  different  shapes  of  openings. 
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Hydrodynamic  Diffusivities  are  linearly  dependent 
on  ambient  current: 
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Figure  67.  Normalized  release  rates  vs.  ambient  current  speed  in  the  range  of  0  <  U  <  0.15  m/s 
(above),  and  the  slopes  for  the  four  different  shapes  of  openings. 
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Figure  68.  Normalized  release  rates  vs.  ambient  current  speed  in  the  range  of  0.15  m/s  <  U  (above), 
and  the  slopes  for  the  four  different  shapes  of  openings. 
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Figure  69.  Normalized  release  rates  vs.  ambient  current  speed  for  rectangle  (4:1  aspect  ratio)  at 
different  orientation  angles. 

4.5.3  Internal  Circulation 

A  FLUENT  simulation  was  conducted  to  obtain  the  internal  eirculation  pattern  inside  the  half 
dome.  The  radius  of  the  cavity,  radius  of  the  breach  hole,  and  the  mean  current  speed  are  R  =  35  mm, 
b  =  2.5  mm,  and  (7=  0.18  m/s,  respectively.  The  steady-state  in-cavity  circulation  patterns  are 
characterized  by  two  jet  steams,  one  pointing  upward  to  the  breach  hole,  and  the  other  which  seems 
to  be  emitted  from  the  hole  and  shooting  downward  (Figure  70).  There  also  appears  to  be  two 
counter-rotating  vortices  on  the  two  sides  of  the  jets.  Similar  circulation  patterns  were  measured  by 
the  PIV  instrument  as  part  of  the  experimental  release  rate  study  (Figure  13),  except  the  downward 
jet  is  not  appreeiable.  The  loeation  and  shape  of  the  upward  jet  and  the  two  eounter-rotating  vortiees 
are  somewhat  different  from  those  seen  in  the  experiment.  Despite  these  differenees,  excellent 
agreement  is  found  between  the  simulation  and  the  experiment  in  terms  of  the  magnitude  of  the 
eireulation  veloeities.  The  variance  of  the  distributed  flow  speed  is  0.012  mm/s  for  simulation  and 
0.014  mm/s  for  the  experiment.  This  is  a  more  important  fact  as  we  argue  that  the  overall  magnitude 
of  circulation  speed  determines  the  in-cavity  mixing  and  thus  the  chemical  release  rate.  It  has  verified 
that  the  in-eavity  eireulation  is  laminar  as  the  overall  eireulation  Reynolds  Number  Rec  =  UcR/  v=  12 
for  the  ease  we  measured  (where  Uc  is  the  average  internal  eireulation  speed,  R  is  the  radius  of  the 
cavity,  and  vis  the  kinematie  viseosity).  Thus,  the  FLUENT  simulation  based  on  the  laminar  model 
is  eonsidered  valid  for  inner  eireulation.  We  ean  also  argue  that  although  the  outer  flow  is  turbulent, 
its  effeet  on  the  release  of  MC  is  confined  to  the  near-breaeh  area.  Given  that  the  size  of  the  breach  is 
small,  the  internal  near-breach  part  is  still  laminar,  i.e.,  in  the  viscous  sublayer  of  the  boundary  layer. 
This  is  partially  verified  by  the  fact  that  the  measured  in-cavity  circulation  is  almost  steady.  These 
two  figures  show  that  the  flow  speeds  between  model  and  measurement  are  of  the  same  order  of 
magnitude,  with  similar  degree  of  variances.  Based  on  both  results,  the  Reynolds  Number  (Ub/u) 

~12,  where  U  is  internal  eirculation  speed  (~5  mm/s),  b  =  2.5  mm,  and  v  is  dynamie  viscosity. 
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Figure  70.  Model-simulated  (above)  and  measured  (below)  internal  circulation. 
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4.5.4  Probabilistic  Modeling 

The  five  variables  governing  the  mass  release  rate  function  (F)  for  a  single  shell,  as  shown  in  Eq. 
(10a)  can  be  prescribed  either  deterministically  or  probabilistically  or  as  a  combination  of  both: 


iTtDC^ 

2D  D  \  \ 

aUb^  juR^  b  R 


(43) 


For  example,  the  sizes  of  the  holes  (with  b  as  the  radius)  on  a  single  shell  will  vary  or  change  in 
time.  The  hole  sizes  may  involve  a  range  or  distribution  function,  which  represents  variability 
associated  with  the  variable  b.  Ambient  current  U  will  also  fluctuate  due  to  the  dynamic  nature  such 
as  tidal  forcing,  wind,  and  freshwater  effects  of  the  marine  environment  in  which  these  shells  reside. 
Therefore,  a  range  or  distribution  function  is  often  used  to  prescribe  a  range  of  expected  currents  in 
the  marine  environment.  Another  source  of  probabilistic  nature  is  the  uncertainty  or  lack  of  data 
associated  with  certain  variables.  Hole  sizes  not  only  involve  variability,  they  often  are  associated 
with  uncertainty  due  to  lack  of  sufficient  data.  To  make  meaningful  predictions,  uncertainty 
associated  with  certain  variables  can  be  prescribed  by  a  probability  function.  Predictions  should  also 
be  interpreted  in  a  probabilistic  fashion  using  an  output  distribution. 

When  the  mass  release  rate  function,  developed  for  a  single  shell,  is  implemented  for  multiple 
shells,  the  number  of  shells  for  each  shell  type  needs  to  be  prescribed,  deterministically  or  probabilis¬ 
tically.  If  it  is  prescribed  probabilistically,  it  can  also  be  processed  in  ways  similar  to  those  for  the 
other  variables.  When  multiple  variables  are  prescribed  probabilistically,  each  with  its  own 
distribution  function,  predictions  would  require  a  large  number  of  calculations  so  that  the  variables 
with  a  probability  distribution  can  be  adequately  sampled.  To  do  this,  we  used  the  Monte  Carlo 
method  to  perform  the  scenario  simulations  (results  described  below).  For  all  of  the  Monte  Carlo 
simulations,  10,000  calculations  were  conducted  for  each  scenario.  Simulation  results  from  the 
10,000  calculations  showed  no  noticeable  difference  from  the  results  of  1,000  calculations. 

4.5.5  Monte  Carlo  Simulation  Results 

The  first  example  scenario  is  for  a  single  24  inches  x  8  inches  round  with  8  kg  of  MC.  Two 
variables  use  probabilistic  input:  the  hole  size  and  ambient  current  speed.  The  hole  radius  is  assumed 
to  be  associated  with  a  normal  distribution  with  a  mean  and  standard  deviation  of  2  and  0.75  cm, 
respectively.  Ambient  current  speed  is  associated  with  a  normal  distribution  with  a  mean  and 
standard  deviation  of  15  and  5  cm/s,  respectively.  Saturation  concentrations  and  dissolution  speed  for 
TNT  were  used  and  hydrodynamic  dispersion  coefficient  was  obtained  from  the  FLUENT  model. 

The  resulting  simulated  release  rates  are  also  associated  with  a  probabilistic  range  (0.01  to  0.033 
mg/s)  and  an  estimated  time  for  the  complete  depletion  of  MC  ranges  from  8  to  12  years  (Figure  71) 

For  the  second  example,  the  hole  size  is  reduced  from  a  2-cm  mean  radius  normal  distribution  to  a 
0.4-cm  mean  radius  log-normal  distribution.  The  predicted  depletion  time  increases  from  10  to  120 
years  (Figure  72). 

The  last  scenario  is  similar  to  the  second  scenario,  except  that  the  number  of  shells  has  a  normal 
distribution  with  a  mean  and  standard  deviation  of  100  and  30,  respectively  Figure  73  [a]  shows  the 
simulated  release  rate  distribution  for  a  single  shell.  When  these  two  numbers  were  multiplied  in  the 
Monte  Carlo  manner,  total  release  rates  were  simulated  and  results  are  shown  in  Figure  73  [c].  Figure 
73  [d]  shows  the  probability  for  each  release  time  or  bin  time  it  takes  to  completely  release  8  kg  of 
MC  from  each  of  100  shells  for  the  scenarios.  And  the  total  cumulative  probability  for  complete 
release  of  all  shells  can  be  obtained  by  integrating  all  the  bins  in  the  curve  (-120  years). 
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Figure  71 .  Monte  Carlo  simulation  for  release  rate  and  time  for  complete  depletion  of  MC;  b  has  a 
normal  distribution. 


Figure  72.  Monte  Carlo  simulations  for  release  rate  and  time  for  complete  depletion  of  MC;  b  has  a 
log-normal  distribution. 
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Figure  73.  Simulated  release  rate  for  group  of  shells,  [a],  prescribed  probability  distribution  for  shell 
number  with  a  mean  of  100  and  standard  deviation  of  30,  [b],  and  simulated  total  release  rate,  [c], 
and  probability  of  release  time,  [d]. 


The  probabilistic  model  is  based  on  the  analytical  release  rate  function.  The  multiple  variables  can 
be  treated  both  deterministically  and/or  probabilistically,  or  a  combination  of  both.  When  the  mass 
release  rate  function  is  implemented  in  the  field,  the  number  of  shells  can  be  added  into  the  function 
as  a  multiplicative  term.  For  real  scenarios,  the  three  variables  that  are  most  likely  associated  with 
probabilistic  estimation,  would  include  number/type  of  shells,  hole  size,  and  ambient  current.  Study 
results  in  this  section  demonstrate  that  the  Monte  Carlo  simulation  is  effective  in  quantifying  the 
effect  on  mass  release  rate  when  multiple  variables  are  associated  with  probability  distributions. 

4.6  EMPIRICAL  STUDY  RESULTS  -  BURIED  SHELL 

Results  from  18  test  runs  described  in  Section  3.4  are  presented  in  Figure  74  through  Figure  91. 
Time  histories  of  the  exit  dye  concentration  Cout  and  log-transformed  Cout  values  are  shown  in  the 
following  figures.  A  log-linear  profile  would  indicate  that  the  model  Eqs.  (28)  and  (29)  proposed  in 
the  above  section  is  valid,  and  that  the  mass  release  velocity  w/can  be  estimated  from  linear 
regression  (note:  sporadic  excursions  or  “spikes”  can  be  found  in  some  datasets,  which  are 
considered  to  be  due  to  sensor  errors). 
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Figure  74.  Fine  sand,  low  seepage  flow  (Ah  =  1”),  crack  hole  diameter  {d^  =  0.21  cm). 
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Figure  75.  Fine  sand,  medium  seepage  flow  (Ah  =  2”),  crack  hole  diameter  (c/b  =  0.21  cm). 


Figure  76.  Fine  sand,  high  seepage  flow  (Ah  =  3”),  crack  hole  diameter  {db  =  0.21  cm). 
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Figure  77.  Fine  sand,  low  seepage  flow  (Ah  =  1”),  crack  hole  diameter  {db  =  0.40  cm). 


Figure  78.  Fine  sand,  medium  seepage  flow  (Ah  =  2”),  crack  hole  diameter  {db  =  0.40  cm). 


Figure  79.  Fine  sand,  high  seepage  flow  (Ah  =  3”),  crack  hole  diameter  {db  =  0.40  cm). 
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Figure  80.  Fine  sand,  low  seepage  flow  (Ah  =  1”),  crack  hole  diameter  {dt,  =  0.56  cm). 


Figure  81.  Fine  sand,  medium  seepage  flow  (Ah  =  2”),  crack  hole  diameter  (db  =  0.56  cm). 


Figure  82.  Fine  sand,  high  seepage  flow  (Ah  =  3”),  crack  hole  diameter  (db  =  0.56  cm). 
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Figure  83.  Coarse  sand,  low  seepage  flow  (Ah  =  1”),  crack  hole  (diameter  =  0.21  cm). 
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Figure  84.  Coarse  sand,  medium  seepage  flow  (Ah  =  2”),  crack  hole  diameter  =  0.21  cm. 
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Figure  86.  Coarse  sand,  low  seepage  flow  (Ah  =  1”),  crack  hole  diameter  db  =  0.40  cm. 


Figure  87.  Coarse  sand,  medium  seepage  flow  (Ah  =  2”),  crack  hole  diameter  db  =  0.40  cm. 
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Figure  90.  Coarse  sand,  medium  seepage  flow  (Ah  =  2”),  crack  hole  diameter  db  =  0.56  cm. 


Figure  91.  Coarse  sand,  high  seepage  flow  (Ah  =  3”),  crack  hole  diameter  db  =  0.56  cm. 
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The  estimated  mass  release  velocity  w/is  plotted  against  the  local  seepage  velocity  Us  in  Figure  92 
for  all  measurements  shown  in  Figure  74  through  Figure  91,  above.  A  linear  relationship  exists 
between  m/ (experimental  mass  release  velocity)  and  (seepage  flow  velocity)  and  the  release 
velocity  has  no  obvious  dependency  on  either  the  crack  hole  diameter  or  on  the  type  of  sand, 
meaning  it  is  not  restricted  by  these  variables.  Linear  regression  (with  intercept  forced  to  zero)  on  all 
of  the  data  indicates  that  Uf=  1.53  Us  with  the  coefficient  of  determination  ^  =  0.96. 


Figure  92.  Relationship  between  the  estimated  mass  release  velocity  {Uf)  and  the  seepage  flow 
velocity  {Us). 

4.6.1  Flux  with  Zero  Seepage  Velocity 

The  limiting  experiment  for  the  case  of  zero  seepage  velocity  was  also  conducted  with  the  same 
container  and  a  breach  hole  diameter  dy,  =  0.56  cm.  A  bucket  (diameter  =  40  cm)  was  filled  with  sand 
and  then  water  to  a  depth  of  20  cm  above  the  sand.  The  breached  container  with  Rhodamine  WT 
(concentration  Co  =  680  ppm)  was  buried  in  sand  at  a  depth  of  ~5  cm.  The  water  above  sand  was 
continuously  and  gently  stirred  (without  causing  seepage  flow  in  the  water)  such  that  once  the  dye 
diffused  from  the  porous  sand  into  the  overlying  water,  it  would  be  adequately  mixed.  Concentration 
of  dye  in  water  was  measured  at  set  time  intervals  (~50  hr)  for  -400  hours  (16  days).  Figure  93 
shows  the  time  series  of  measured  dye  concentrations  in  water. 

The  remarkably  linear  variation  of  concentration  with  time  indicates  a  quasi-constant  release  rate 
of  dye  from  the  container  over  the  2-week  period.  For  the  release  rate  model,  F  =  CoAuf,  the 
estimated  mass  release  velocity  in  this  case  is  Uf=  0.00005  (cm/s),  -2  orders  of  magnitude  smaller 
than  the  case  with  seepage  flow  exhibiting  the  lowest  mass  release  velocity  of  all  cases  tested  (m^  = 
0.023  mm/s). 
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Figure  93.  Change  of  concentration  in  water  above  the  sand. 

4.7  MODELING  STUDY  RESULTS  -  BURIED  SHELL 

Figure  94  compares  the  Darcy  velocity  calculated  by  the  model  to  the  measured  velocity  for  the 
fine  and  coarse  sand  under  the  three  different  head  gradient  scenarios.  Model-simulated  pore-water 
velocities  mimicked  measured  values  for  all  of  these  scenarios. 


Figure  94.  Comparison  of  simulated  and  measured  Darcy  velocities  at  three  head  pressure 
differences  for  fine  and  coarse  sands. 
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4.7.1  Transport  of  TNT  in  Sediment 

As  discussed  previously,  the  release  rate  of  MC  from  a  single  breached  shell  (e.g.,  Eq.  (10a)  or  Eq. 
(43))  depend  on  five  parameters,  which  includes  ambient  current  (U),  hole  size  (b),  cavity  volume 
(R),  dissolution  rate  constant  (p),  and  mixing  coefficient  (D).  These  can  be  grouped  into  two 
Reynolds  numbers,  namely  Ub/D  and  pR/D.  We  have  demonstrated  that  the  MC  release  rate  is 
limited  by  the  combination  of  these  two  Reynolds  numbers; 

(i)  Current-Controlled  Release  (pR  /D  »  Ub/D) 

When  the  dissolution  speed  (p)  is  much  larger  than  the  ambient  current  (U)  or  the  cavity  (R)  is 
much  larger  than  the  hole  (b),  MC  solutions  inside  the  cavity  are  near  saturation  and  the  release  rate 
of  MC  is  controlled  by  the  ambient  current  speed  U  and  the  hole  size. 

(ii)  Dissolution  Rate-Controlled  Release  (pR  /D  «  Ub/D) 

For  this  case,  release  of  MC  is  controlled  by  the  dissolution  rate  (p)  since  ambient  current  and/or 
the  hole  is  so  large  that  MC  solution  dissolved  from  the  solid  phase  is  dispersed  out  of  the  breached 
shell  immediately. 

For  both  cases,  the  mixing  coefficient,  D,  plays  a  significant  role  in  the  release  rate.  We  previously 
demonstrated  that  for  the  water  column,  current  speed  is  in  general  larger  than  0.1  cm/sec,  which 
according  to  the  mixing  coefficient,  developed  empirically  can  be  estimated  as 

—  =  omiub, 

P 


where  |3  is  an  empirical  constant. 

Generally,  ambient  current  and  its  induced  mixing  coefficient  is  the  primary  driving  force  for  the 
release  rate.  For  release  in  sediment,  MC  is  subject  to  mixing  by  low  advective  flow  and  molecular 
diffusivity.  As  the  advection  of  the  pore-water  decreases,  the  current-induced  mixing  would  be 
reduced  accordingly  to  a  level  comparable  to  molecular  diffusivity,  as  shown  in  Figure  95.  As  an 
extreme  condition,  when  advection  of  pore  water  becomes  zero  (stagnant  flow),  fate  and  transport  of 
MC  in  sediment  is  controlled  by  molecular  diffusivity  alone. 

For  example,  in  the  case  of  TNT,  molecular  diffusivity  becomes  important  when  ambient  current 
speed  is  reduced  below  ~10'^  cm/s  (Figure  95).  This  low  current  speed  can  be  considered  to  be 
stagnant  or  below  detection  limits  for  surface  water.  However,  for  marine  sediment,  this  low  pore 
water  flow  is  within  the  range  for  pore  water  advection.  For  sediment,  molecular  diffusivity  needs  to 
be  corrected  for  porosity  (or  tortuosity)  of  the  sediment  medium; 


Dfiff  -  Dmoi  sot,  (44) 

where  a  is  the  porosity  and  a  is  an  empirical  constant  (-1.2  to  1 .8)  validated  by  the  empirical  study, 
and  Dmoi  is  the  molecular  diffusivity  of  TNT  (6.54  x  10'®  cm^/s)  (Lynch,  Brannon,  and  Delfino, 
2002a).  Using  the  diffusivity  and  pore-water  flow  with  speed  of  0.00001,  0.0001,  and  0.001  cm/s, 
respectively,  the  release  rate  and  TNT  plume  are  modeled  under  three  scenarios  with  different 
relative  magnitudes  between  the  advection  and  diffusivity  (molecular)  extremes. 
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Advection-Diffusivity  Transport  of  MC  in  Sediment  and  Water 
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Figure  95.  Range  of  advection  and  diffusivity  {e.g.,  molecular)  in  water  and  sediment  (pore  water). 

4.7.2  Advection-Diffusivity  Transport  of  MC  in  Sediment 

Transport  of  MC  in  sediment  is  governed  by  the  advection-dispersion  transport  equation,  which  is 
the  same  fundamental  equation  as  used  for  transport  in  the  water  column.  In  the  latter  case,  advection 
by  ambient  currents  dominate  the  mixing  and  transport  of  MC.  As  discussed  above,  for  transport  in 
sediment,  both  advection  and  diffusivity  (molecular)  need  to  be  considered,  in  particular,  for  low  or 
stagnant- water  flows. 

Figure  96  through  Figure  98  show  simulated  dispersive  behaviors  of  a  substance  (e.g.,  MC), 
located  at  (0,0)  with  concentration  of  500  units  initially  (t  >  0)  at  time  =  23  hours,  7th  day  and  30th 
day,  respectively.  A  horizontal  diffusivity  Dx  =  Dy  =  1.18X10-6  cm2/s  for  the  x-y  plane  and  a 
vertical  diffusivity  =  0.5  Dx(y)  for  the  z-direction  are  used.  Concentrations  in  the  three  vertical 
planes,  z  =  -5,  0  and  5  cm,  respectively,  are  shown.  The  figure  shows  that  on  the  z  =  0  plane, 
diffusion  of  concentration  at  the  10%  of  the  boundary  condition  (500  units  at  (0,0))  expand  to  a 
radius  of  ~5  cm  and  ~1  cm  at  the  end  of  7th  day  and  30th  day,  respectively.  At  both  z  =  -5  cm  and  5 
cm  planes,  no  appreciable  concentrations  are  reached  at  t  =  23  hours.  These  different  concentrations 
are  within  the  same  order  of  magnitude  as  those  measured  in  the  laboratory,  discussed  previously. 

Figure  99  shows  the  advection-diffusion  transport  for  different  relative  magnitudes  of  advection 
and  diffusivity.  The  three  pore-water  seepage  velocities  are  0.001,  0.0001,  and  0.00001  cm/s,  and  the 
molecular  diffusivities  in  the  x-y  and  z  planes  are  the  same  as  described  above  (Figure  96  through 
Figure  98).  It  is  hypothesized  that  the  three  seepage  velocities  are  selected  to  cover  a  flow  regime  for 
advection  dominated  flow  (U  =  0.001  cm/s)  and  a  flow  regime  governed  equally  by  advection  and 
diffusivity  (U  =  0.00001  cm/s).  Flow  advection  dictates  the  length  of  the  plume,  whereas  diffusivity 
dictates  the  expansion  (width)  of  the  plume.  In  Figure  99,  the  characteristics  of  advection  and 
diffusivity  effects  on  the  transport  can  be  visualized  clearly  in  the  model-predicted  contours  of  the 
plumes. 
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Figure  96.  Dispersion  of  TNT  in  sediment  by  molecular  diffusivity  (Dx  =  Dy  =  1.18  x  10"®  cm^/sec,  Dz 
0.5  Dx(y))  at  time  =  23  hours  from  initial  condition  (continuous  source  for  t  >  0,  with  concentration  of 
500  units  at  (0,0)). 
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Figure  97.  Dispersion  of  TNT  in  sediment  by  molecular  diffusivity  (Dx  =  Dy  =  1 .18  x  10'®  cm^/sec, 

Dz  =  0.5  Dx(y))  at  time  =  7th  day  from  initial  condition  (continuous  source  for  t  >  0,  with  concentration 
of  500  units  at  (0,0)). 
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Figure  98.  Dispersion  of  TNT  in  sediment  by  molecular  diffusivity  (Dx  =  Dy  =  1.18  xIO'®  cm^/sec, 

Dz  =  0.5  Dx(y))  at  time  =  30th  day  from  initial  condition  (continuous  source  for  t  >  0,  with  concentration 
of  500  units  at  (0,0)). 
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Figure  99.  Transport  of  TNT  in  sediment  by  molecular  diffusivity  (D  =  1.18  x  10'®  cm^/s)  and 
advection  (U  =  0.001,  0.0001  and  0.00001  cm/s,  for  top,  middle,  and  bottom  figures,  respectively) 
from  initial  condition  (continuous  source  fort  >  0,  with  concentration  of  100,000  units  at  (0,0)) 
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5.  CONCLUSIONS  AND  IMPLICATIONS  FOR  FUTURE 
RESEARCH/IMPLEMENTATION 

Theoretical  and  numerical  studies  were  conducted  for  the  release  of  munitions  compounds  from  a 
single  breached  shell  under  various  hydrodynamic  and  shell  integrity  conditions.  Analytical  solutions 
were  derived  for  the  release  rate  based  on  the  assumptions  that  ambient  currents  are  tangential  over 
the  shell  surface,  and  that  the  hole  is  small  in  size  compared  to  the  cavity  inside  the  shell.  These 
theoretical  analyses  and  numerical  simulations  have  demonstrated  that  the  flux  of  MC  solution  can 
be  expressed  explicitly  as  a  function  of  five  primary  variables;  current  speed  U,  hole  area^,  saturated 
concentration  Cs ,  dissolution  speed  (p),  and  hydrodynamic  dispersion  coefficient  D  (Eq.  (10a)).  The 
analytical  release  rate  solutions  were  validated  by  the  model  results  using  FLUENT. 

5.1  RELEASE  RATE 

The  validated  results  show  that,  in  a  dynamic  estuarine  water  body,  release  rates  can  experience 
processes  that  cover  wide  regimes  including  two  extreme  limiting  conditions.  For  the  first  condition, 
ub^  «  pR^,  the  release  rate  is  limited  by  the  ambient  current  speed.  MC  dissolution  in  the  shell 
cavity  occurs  at  a  faster  rate  than  the  ability  of  MC  to  escape  through  the  hole  as  a  function  of  the 
ambient  current  (hydrodynamic  dispersion).  As  a  result,  the  solution  in  the  cavity  is  always  near 
saturation.  For  the  condition,  ub^  »  pR^,  where  the  release  rate  is  limited  by  dissolution  speed,  MC 
dissolution  in  the  cavity  is  much  slower  than  the  transport  of  MC  solution  through  the  hole.  For  this 
scenario,  the  solution  in  the  cavity  is  not  near  saturation  and  the  retention  time  of  the  MC  solution 
inside  the  breached  shell  (cavity)  is  small,  resulting  in  more  mass  per  unit  time  escaping  from  the 
shell. 

When  the  thickness  of  the  shell  (/?)  is  considered,  the  general  release  rate  retains  the  same 
functional  dependence  on  the  five  variables  above.  However,  the  shell  thickness  imposes  an  effect  on 
the  release  rate,  which  can  be  expressed  explicitly  as 

F  =  kU  ,  (45) 

where  /c  is  a  coefficient  which  assimilates  the  effect  of  the  geometry  of  the  breach  and  shell 
thickness.  If  the  breach  has  a  circular  shape,  the  value  of  k  is  between  0.015  and  0.035  according 
to  numerical  simulations,  in  comparison  with  value  of  0.023  when  the  thickness  is  ignored.  This 
empirical  model  suggests  that  the  MC  release  rate  grows  linearly  with  U,  A,  and  Cs,  while  it  is 
independent  of  the  cavity  size  R.  This  model  is  only  valid  if  the  breach  size  is  relatively  small 
compared  to  the  cavity  inside  the  shell  (b  «  R)  and  the  shell  thickness  is  large  compared  to  the 
hole  radius  (h>  b). 

A  series  of  scenarios  were  designed  and  conducted  in  Phase  III  to  study  the  release  rate  of  MC 
from  breached  shells  in  marine  environment.  Based  on  the  semi-analytical  release  rate  function 
developed  in  the  Phase  I  study,  extensive  experiments  and  corresponding  modeling  studies  were 
conducted  to  examine  the  mass  release  rate  of  a  surrogate  dye  from  a  semi-spherical  cavity  through  a 
circular  hole.  The  experiments  investigated  the  effects  of  breach  hole  sizes,  shapes,  and  orientations 
on  transport  out  of  the  breached  shell  into  the  external  water.  Of  the  five  variables  that  define  the 
release  rate,  the  hole  size  and  the  ambient  current  speed  appear  to  be  the  two  most  important 
parameters.  As  such,  the  internal  circulation  magnitudes  were  also  measured,  analyzed,  and 
characterized.  The  mass  release  rate  is  considered  to  be  proportional  to  the  local  shear  velocity,  the 
mass  concentration  in  cavity  and  the  area  of  the  breach.  Dimensional  analysis  indicates  that  the  mass 
release  rate  can  be  defined  by  a  characteristic  speed  Uf,  which  is  defined  as  the  mass  release  speed. 
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i.e.,  F  =  UfACo,  where  A  is  the  area  of  the  crack,  and  Co  is  the  concentration  of  mass  within  the  cavity. 
The  mass  release  speed  would  then  depend  on  the  shear  velocity  above  the  crack,  u*,  the  size  of  the 
crack  and  the  hole,  as  well  as  the  kinematic  viscosity  of  water. 

Experimental  study  in  this  report  entailed  the  variables  with  values  with  the  following  ranges:  Uf  = 
0.4-2  u*  for  b  =  0.5  -  25  mm  and  F  =  1  ~  5  cm.  It  also  shows  that  there  is  no  clear  evidence  of  a 
relationship  between  the  ratio  of  Uf/u*  and  the  hole  radius  b,  or  the  Reynolds  Number  u*b/v. 

However,  Uflu*  is  largely  determined  by  the  cavity  size  R,  and,  thus,  the  Reynolds  Number  u*R/v.  An 
extrapolation  of  the  experimental  results  indicates  that  Wf/w*  is  approximately  3  to  4  for  R  =  10  cm, 
which  is  the  case  for  our  previous  FLUENT  simulation.  Numerical  simulation  had  shown  that  Wf/w*  in 
that  case  was  about  1 . 

A  linear  model  between  m/m*  and  the  cavity-based  Reynolds  number  Rcr  seems  to  be  plausible 
according  to  this  series  of  experiments.  The  strong  dependency  on  the  cavity  size  R  was  not  explicitly 
considered  in  our  simple  theoretical  model,  that  is. 


F  =  UC{b)nb^ 


2D  Dll' 

aUb^  aR^  b  R 


(46) 


However,  as  we  assume  infinite  dissolving  speed  and  large  Rib  ratio,  the  mass  release  rate  can  be 
described  as 


2D  ]_  ■  (47) 

aUb^  b 

Additionally,  the  dependency  on  R  could  be  linked  through  the  in-cavity  mixing  coefficient  D. 

TNT  mineralization  occurs  by  multiple  pathways,  including  transformation,  bio-degradation  and 
mineralization.  Mineralization  rates  of  TNT  as  reported  are  based  on  laboratory  studies,  including 
those  of  Dr.  Mike  Montgomery  of  NRL  (Montgomery  et  al,  2008).  The  mineralization  rates, 
conducted  and  measured  in  the  laboratory  experiments  appear  to  be  high  and  do  not  follow  zero-  or 
first-order  criteria.  With  mineralization  rates  this  high,  simulated  TNT  concentrations  should  be  very 
low.  In  this  study,  we  only  consider  a  reported  mineralization  process  by  mineralization.  Thus,  the 
reported  1.5%  TNT  mineralization  rate  over  a  24-hour  period  assuming  concentrations  of  TNT  of  0 
to  400  ng/L  (water)  or  ng/Kg  (sediment).  Assuming  a  first-order  mineralization  (an  unverified 
assumption  and  extrapolation  from  laboratory  data),  this  translates  into  in  a  half-life  of  0.0151/day. 

5.2  MC  TRANSPORT 

Transport  processes  for  TNT  were  simulated  using  both  a  fine-scaled  hydrodynamic  and  transport 
model,  TRIM2D,  and  a  box-scaled  model,  WASP7.  Both  of  these  models  provide  meaningful  results 
that  are  mutually  consistent  and  complementary,  demonstrating  their  utility  for  developing  an 
understanding  of  munitions  constituent  transport  in  marine  environment. 

For  the  conservative  scenario,  where  TNT  release  loads  enter  the  inner  San  Diego  Bay  (Box  #1  for 
WASP)  water  column  with  no  mineralization  or  exchange  of  mass  with  the  sediment,  results  from 
both  models  compare  well,  demonstrating  the  consistency  of  both  models.  Results  for  this  scenario 
provide  conservative  estimates  of  TNT  in  San  Diego  Bay. 
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With  the  use  of  reported  soil  partitioning  coefficients  for  TNT,  which  are  low,  only  a  small 
fraction  of  TNT  mass  is  distributed  into  the  particulate/sediment  phase.  WASP-simulated  TNT 
concentrations  in  sediment  (<  0. 1  ng/L  for  San  Diego  Bay  and  <  1  ng/L  for  Elizabeth  River)  are 
much  lower  in  this  scenario  than  the  simulated  TNT  concentrations  in  the  water  column  (<  45  ng/L). 
Therefore,  the  particulate  phase-associated  TNT  in  both  water  column  and  sediment  should  be  low. 
These  model  results  are  based  on  the  assumption  of  equilibrium  partitioning  between  particulate  and 
aqueous  phase  of  TNT,  which  is  considered  a  reversible  process.  However,  as  discussed  by 
Brannon,  it  was  observed  that  TNT  is  lost  from  the  water  column  at  relatively  fast  rates  and  sediment 
is  a  sink  for  TNT  in  the  water  column.  This  can  be  attributed  to  continuing  removal  of  the  compound 
from  the  water  by  degradation,  transformation,  or  covalent  bonding  to  sediment  components  such  as 
organic  carbon  (Brannon  et  ah,  2005).  Mostly  likely,  these  processes  are  not  reversible,  contrary  to  a 
partitioning  process,  which  is  considered  reversible. 

5.3  SETTLING  AND  RESUSPENSION 

When  the  exchange  of  TNT  mass  between  water  column  and  sediment  is  considered,  WASP7  is 
used  to  account  for  these  mass-exchange  processes,  including  partitioning,  settling,  and  resuspension, 
in  additional  to  the  hydrodynamic  transport.  Partitioning  of  TNT  onto  TOC  with  TSS  as  the  transport 
vehicle  describes  the  particulate  phase-associated  behavior  of  TNT.  Mass  exchange  of  TNT  across 
the  water/sediment  interface  takes  place  through  settling  and  resuspension  of  TSS.  By  settling,  the 
particulate  phase  of  TNT  moves  into  the  sediment  bed,  which  is  assumed  to  be  homogenous  over  the 
top  mixing-layer  thickness.  Sediments,  and  the  particulate  phase  of  TNT,  are  considered  mobile  and 
are  transported  back  into  the  water  column  through  the  resuspension  processes.  However,  the 
empirical  study  conducted  by  Pennington  indicated  that  sorption  of  TNT  to  TOC  in  water  could  not 
reach  a  steady-state  condition,  implying  that  sorption  of  TNT  may  not  follow  a  reversible  process, 
which  is  the  underlying  assumption  for  use  of  partitioning  coefficient.  Accordingly,  for  this  study, 
partitioning  coefficients  reported  for  soils  were  used  to  simulate  TNT  in  both  the  water  column  and 
sediment. 

Generally,  literature  reports  describing  gross  settling  and  resuspension  rates  (velocities)  in  rivers 
and  harbors  are  not  adequate.  However,  field  data,  such  as  net  erosion  or  deposition  over  extended 
time,  can  be  used  to  deduce  the  order  of  magnitude  for  each  individual  process.  In  this  study,  we  used 
Stokes  Law  to  estimate  settling  velocity  and  the  resuspension  velocity  is  calculated  based  on  the 
assumption  that  there  is  no  net  erosion  or  deposition  throughout  San  Diego  Bay.  This  baseline 
assumption/approach  is  taken  for  simplification  purposes  in  the  modeling  study.  As  more  field  data 
become  available,  estimates  with  higher  certainty  can  be  made  for  the  settling  and  resuspension 
velocities. 

For  the  release  rate  of  8  kg/year  from  a  single-shell,  it  is  estimated  that  TNT  concentrations  in  the 
water  column  are  ~40  ng/L  for  both  estuaries.  Simulated  TNT  concentrations  in  sediment  for  San 
Diego  Bay  are  ~0.1  ng/L,  and  ~1  ng/L  for  the  Elizabeth  River,  a  one  order  of  magnitude  difference. 
This  is  because  the  two  systems  have  significantly  different  TSS  characteristics  and  behaviors  (Table 
19). 

For  example,  TSS  concentrations  in  San  Diego  Bay  and  Elizabeth  River  are  2  mg/L  and  15  mg/L, 
respectively,  a  factor  of  7X.  Settling  velocities  are  12  m/day  for  San  Diego  Bay  and  50  m/day  for 
Elizabeth  River,  a  difference  of  a  factor  of  4X.  Resuspension  velocity  for  San  Diego  Bay  is  2.6  x  10'^ 
m/day,  compared  to  Elizabeth  River  at  3.0  x  10"^  m/day,  a  difference  of  three  orders  of  magnitude. 
Therefore,  the  exchange  of  TSS  and  associated  TNT  between  the  water  column  and  sediment,  via 
settling,  resuspension,  and  partitioning,  is  much  more  active  in  Elizabeth  River  than  in  San  Diego 
Bay,  resulting  in  higher  predicted  TNT  concentrations  in  sediment. 
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Table  18.  Comparison  of  TSS,  settling  velocity,  and  resuspension  velocity  for  San  Diego  Bay  and 
Elizabeth  River. 


San  Diego  Bay 

Elizabeth  River 

Settling  Velocity  (meters/day) 

12 

50 

Resuspension  Velocity  (meters/day) 

2.6x10'^ 

3.0x10'' 

Total  Suspended  Solid  (mg/L) 

2 

15 

Transport  and  flushing  of  munition  compounds  in  estuaries  is  primarily  governed  by  hydrody- 
namies,  ineluding  tides  and  freshwater  inflows.  We  have  examined  the  transport  and  flushing 
charaeteristics  for  TNT  in  two  estuarine  seenarios:  San  Diego  Bay,  whieh  has  only  tidal  flushing,  and 
Elizabeth  River,  which  has  both  tidal  and  riverine  flushing.  These  examples,  chosen  for  their 
respectively  unique  estuarine  characteristics,  should  represent  most  transport  and  flushing  behaviors 
encountered  in  DoD  harbors.  . 

5.4  SCALING  AND  INTERPRETATION 

An  energy  balanee  analysis  was  undertaken  based  on  the  sketeh  shown  in  Figure  100.  The  overall 
magnitude  of  in-cavity  circulation  velocity  is  denoted  as  Mm-  The  eharacteristie  veloeity  near  the 
breach  can  be  considered  as  the  shear  veloeity  u*.  Taking  the  cavity  as  the  control  volume,  the  energy 
influx  into  the  eavity  should  be  in  balanee  with  energy  dissipation  in  the  volume  under  an 
equilibrium  state.  The  energy  influx  should  scale  with  u*^b^  and  the  dissipation  should  scale  with 
{umlR)l^  =  UmR^.  If  this  scaling  is  right,  the  mixing  coefficient  can  be  scaled  as 

(48) 

thus  D  becomes  relevant  in  Equations  (47)  and  (47). 

Experimental  investigation  on  the  sealing  between  Mm  and  u*  is  underway,  whieh  may  explain  the 
observed  dependeney  of  Mf  on  the  cavity  size  R. 


Figure  100.  A  sketch  of  an  energy  balance  model  to  scale  the  in-cavity  circulation  to  external  shear 
velocity. 

Overall,  assuming  equal  opening  areas,  different  shapes  of  the  breach  hole  will  generate  different 
mass  release  rates,  which  are,  however,  small  to  moderate  (Figure  66,  Figure  67,  and  Figure  68).  The 
size  of  the  hole  (area)  has  more  effeet  on  the  mass  release  rate  than  the  shape  of  the  hole.  The  same 
eonclusion  is  made  for  the  orientation  of  eracks,  which  are  simulated  as  two  rectangles  with  the 


104 


aspect  ratios  of  1:4  and  1:10,  respectively  (Figure  69).  These  findings  indicate  that  hydrodynamic 
mixing  (i.e.,  release  velocity)  is  governed  more  by  the  magnitude  of  the  ambient  current  than  by  the 
direction  of  the  current  relative  to  the  breach  holes.  In  reality,  formation  of  the  breach  holes,  and 
cracks  by  corrosion  can  be  a  random  process  in  terms  of  their  shapes  and  orientations,  which  is  not  as 
important  for  the  mass  release  rates  as  is  the  total  area  of  the  breach  hole. 

While  the  mass  release  rate  for  a  single  shell  with  a  breach  hole  has  been  studied,  there  are 
inherent  uncertainties  for  estimation  for  release  from  more  than  one  shell.  The  semi-empirical  release 
function,  which  has  been  validated  by  both  the  empirical  and  modeling  studies,  can  be  used  to 
estimate  the  release  rates  of  MC  (surrogate  dye  in  the  study)  from  a  breached  shell  under  various 
cases  of  shell  integrity  and  ambient  hydrodynamic  conditions.  Uncertainty  for  the  release  rate 
includes  two  sources.  The  first  sources  of  uncertainty  are  variabilities  associated  with  the  five 
variables  in  the  release  rate  function.  For  example,  a  group  of  breached  shells  may  have  breach  holes 
with  a  variety  of  various  sizes.  The  ambient  hydrodynamic  current  speed  may  also  fluctuate  with 
different  tidal  and  flow  conditions  under  which  the  breach  shells  reside.  As  such,  probabilistic 
modeling  simulations  were  conducted  (Section  0).  In  these  simulations,  the  size  of  the  breach  hole 
and  the  ambient  current  are  allowed  to  possess  probability  distributions  due  to  variabilities  associated 
with  them.  The  Monte  Carlo  method  was  used  to  quantify  the  effects  of  these  two  variables  on  the 
mass  release  rate  in  a  probabilistic  manner  (Figure  71,  Figure  72,  and  Figure  73).  The  same 
probabilistic  modeling  method  can  be  used  to  study  the  uncertainties  associated  with  lack  or 
inadequacy  of  the  knowledge  or  data  associated  with  the  variables  for  the  mass  release  rate. 

5.5  SIMULATING  BREACHED  SHELL  SCENARIOS 

Initially,  MC  release  into  the  water  column  from  breached  shells  in  marine  environments  was 
studied  using  three  complimentary  approaches:  analytical,  empirical,  and  numerical  study  compo¬ 
nents.  The  study  was  based  on  the  conceptual  model  shown  in  Figure  1  with  MC  release  from  a 
breached  shell.  Two  initial  nominal  scenarios  were  studied  and  results  validated,  including  a 
breached  shell  sitting  on  top  of  the  sediment  with  hydrodynamic  actions  imposed  by  ambient  current, 
and  a  low  order  detonation  (MC  fully  exposed  to  ambient  current).  We  have  also  addressed  a  third 
nominal  case,  a  breached  shell  buried  in  sediment,  for  which  a  study  of  MC  release  from  breached 
shells  buried  in  sand  has  been  conducted  and  results  reported  herein.  In  a  manner  similar  to  the  water 
column  release  studies,  fluorescent  dye  experiments  were  conducted  in  a  laboratory  sediment  channel 
contained  within  a  constructed  sand  flume  box.  Two  types  of  sand  were  used  as  the  porous  medium 
to  simulate  the  sediment  bed  of  a  coastal  ocean.  The  mean  grain  diameter  of  the  two  sands  were  ds  = 
150  pm  and  500  pm.  A  small  bottle  of  fluorescent  dye  solution  with  a  breach  hole  in  the  cap  (hole 
diameter  =  0.21,  0.40,  and  0.56  cm)  was  buried  in  the  sediment  channel  and  the  dye  was 
transported  by  seepage  flow  under  a  maintained  pressure  gradient.  Eighteen  experiments  were 
conducted  to  simulate  shells  buried  in  sediment,  which  included  two  sand  sizes  (fine  and  coarse 
sand),  three  pore-water  flow  regimes  (low,  medium,  and  high),  and  three  hole  sizes  (0.2,  0.4,  and 
0.56  cm). 

Measured  pore-water  flow  rate  (speed)  compared  well  with  the  model-simulated  values  for  all  of 
these  cases,  confirming  that  results  from  the  empirical  and  numerical  models  are  consistent,  and 
should  thus  be  accurate.  Release  rates  were  determined  by  measuring  the  release  dye  concentrations 
at  the  outlet  of  the  flume  box  over  time.  The  measured  time  series  of  the  dye  concentration  compared 
well  with  the  analytical  solution  for  all  18  cases  and  follows  a  first-order  decay  rate.  The  good 
comparison  of  the  pore-water  advection  over  the  low-,  medium-,  and  high- flow  regimes,  coupled 
with  the  dye  release  rate  provides  evidence  that  processes  controlling  MC  release  from  breached 
shells  buried  in  sand  have  been  successfully  parameterized,  simulated,  and  validated.  A  first-order 
depletion  model  can  be  used  to  explain  the  observed  time  history  of  the  exit  dye  concentration.  This 
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model  is  then  applied  to  validate  the  mass  release  rate  model, 

F  =  CqAuj  , 

where  Co  is  the  mass  concentration  near  the  breach,  A  is  the  area  of  the  breach,  and  w/is  a  mass 
release  velocity.  Experimental  results  indicate  that  w/is  linearly  proportional  to  the  seepage  velocity 
near  the  breach  hole  and  1.53  times  the  seepage  velocity  (m,).  The  study  of  mass  release  rate  from 
breach  holes  subject  to  water  flow  shows  that  Uf=  0.4~2  u*,  where  u*  is  the  shear  velocity  of  the 
mean  flow.  In  both  water  column  flow  and  sediment  porous  seepage  flow,  mass  release  velocity,  Uf, 
is  linearly  proportional  to  the  index  velocity,  which  is  the  friction  velocity,  u*,  and  the  seepage  flow 
velocity,  Us,  for  water  column  flow  and  sediment  seepage  flow,  respectively.  It  is  important  to 
observe  that  the  proportionality  constants  of  1.53  (Figure  92)  for  sediment  seepage  flow  and  0.4to  2.0 
for  water  flow  (Wang  et  al,  2007)  are  comparable  in  magnitude. 

5.6  FUTURE  RESEARCH/IMPLEMENTATION 

Much  understanding  has  been  obtained  about  the  fate  and  transport  of  MC  leached  from  breached 
shells  in  both  the  water  column  and  sediment.  When  the  shell  is  buried  in  sediment,  fate  and  transport 
of  MC  released  from  the  shell  would  be  slowed,  occurring  on  a  time  scale  much  longer  than  the 
modeled  transport  in  the  water  column.  This  slow,  long-term  process  needs  to  be  considered  in 
conjunction  with  the  possibility  of  accumulated  effects  from  multiple  buried  shells.  A  quantitative, 
model-based  treatment  of  these  factors  will  be  required  to  understand  the  significance  of  fate  and 
transport  for  these  scenarios,  which  may  be  reflective  of  a  more  realistic  or  real-world  scenario. 
Future  experimental  and  monitoring  studies  are  needed  to  further  understand  the  behavior  of  TNT  in 
marine  waters,  in  particular  the  ability  to  evaluate  dynamic  exchange  and  transformations  of  TNT 
and  degradation  products  specific  to  known  decay  pathways  appear  to  be  limited.  This  type  of  effort 
should  include  a  systematic  sensitivity/uncertainty  evaluation  and  selective  experimental  and 
monitoring  validation  of  the  developed  critical  model  components. 

Finally,  MC  transport  modeling  through  sediment  surrounding  a  buried,  breached  shell  may 
include  the  presence  of  naturally  occurring  organic  matter  (NOM).  Such  organic  matter  may  impact 
both  MC  transport  and  related  degradation.  While  the  study  of  NOM  as  a  factor  in  attenuating  MC 
fate  and  transfer  was  not  included  in  this  effort,  this  issue  can  be  considered  by  the  shell  model  in 
two  different  ways;  (a)  as  dissolved  NOM  associated  chemically  with  MCs  (as  MC-NOM),  and  (b) 
as  MC  partitioning  onto  particulates  of  NOM.  The  shell  model  can  treat  MC-NOM  in  a  similar 
fashion  as  that  for  MC  in  case  (a),  and  it  is  readily  implementable  to  simulate  dissolved  MC-NOM  in 
a  way  similar  to  that  for  MC  in  case  (b).  However,  MC-NOM  transport  and  MC  transport  are  similar 
enough  that  explicit  treatment  may  not  be  needed.  The  MC  partitioning  (b)  can  also  be  treated  as  part 
of  sedimentation  and  resuspension,  a  capability  which  is  included  in  the  model. 
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